Modeling of solid-supercritical fluid 

EQUILIBRIA USING CUBIC EQUATIONS OF STATE: A 

UNIFIED APPROACH 


A Thesis Submitted 

in Partial Fulfillment of the Requirements 
for the Degree of 
Master of Technology 


by 

B. Chender Rao 


to the 

Department of Chemical Engineering 
Indian Institute of Technology, Kanpur 


February 1996. 






• CERTIFICATE 

It is certified that the work contained in the thesis titled Modeling of Solid- 
Supercritical Fluid Equilibria Using Cubic Equations of State: A Unified 
Approach, by B. Chender Rao , has been carried out under my supervision and that 
this work has not been submitted elsewhere for a degree. 



Dr. R.P. Singh 


Professor 

Department of Chemical Engg. 
IIT Kanpur 


22 February, 1996 



Acknowledgement 


I express my deepest sense of gratitude and appreciation for my thesis supervisor 
Dr. R.P. Singh for his guidance, encouragement, criticisms and advice made during 
this work. 

It is a pleasure to thank my friends and classmates who made the stay at IIT Kanpur 
a memorable one. 

My special thanks to the faculties and authorities of IIT Kanpur who make this in- 
stitution an excellent learning place. 

I take this opportunity to thank Engineers India Ltd., for the financial support pro- 
vided by them during this work. 


B. Chcnder Rao 



Contents 


Abstract iv 

List of Figures vi 

List of Tables vii 

Nomenclature viii 

1 Introduction 1 

2 Literature Survey 4 

3 Thermodynamic Modeling For Solid-Fluid Equilibria 10 

3.1 Equation of State Modeling of Solubilities 10 

3.2 Cubic Equations-of-State 11 

3.3 Fugacity Coefficient from Cubic EOS 12 

3.4 A Simplified Cubic EOS Model for Solubility 14 

4 Results and Discussion 17 

4.1 Physical Properties 17 


4.2 Data Reduction and Parameter Estimation Criteria 


19 



CONTENTS iii 

t 

4.3 Results and Discussion 20 

4.4 Conclusions 61 

Bibliography 62 

Appendix 69 

A Program Listing of the Rigorous Model 70 

B Program Listing of the Simplified Model 93 




Abstract 


The appropriate mathematical modeling of the thermodynamic behaviour of the 
solute-solvent system is important for effective design of the extraction and sepa- 
ration units of a supercritical fluid extraction (SFE) process. In the present work, 
several combinations of cubic equations of state (EOS) and mixing rules have been 
studied, computing binary interaction parameters (BIPs) and absolute average errors 
in solubility for binary systems of 14 solid solutes and 5 supercritical solvents at dif- 
ferent temperatures . Equations of state chosen for this study were Redlich-Kwong 
(R.K), Soave (SRK) and Peng- Robinson (PR), whereas the mixing rules chosen were 
quadratic and cube-root using a single BIP, k x2 or two BIPs, k x2 and l x2 . As expected, 
the introduction of the second BIP, l i2 , greatly improves the accuracy of the calcu- 
lated solubility values. The results are highly satisfactory, considering the uncertainity 
in experimental measurements of low solubility values as well as the uncertainity of 
the T,., P r and u values of the solute. The comparison of overall percent, average 
errors indicate that, globally, no distinction can be made either among the EOS or 
the mixing rules. However, differences among the EOS are found when the errors for 
individual systems are considered. Promising results have been obtained also using 
a recently proposed simplified model with each of the three cubic EOS studied. In 
this method the infinite-dilution fugacity coefficient of the solute in the supercriti- 
cal fluid phase is computed using the infinite-dilution compressibility factor, and the 
EOS parameters a and ft, defined by a = \Ja 2 /a x and /? = b 2 /b x , are treated as 
adjustable parameters; the conventional iterative procedure is thus substituted by a 
once-through calculation. Optimum solute- to-solvent parameter ratios a and p and 
percent average errors have been computed for binary systems of 15 solid solutes and 
5 solvents at different temperatures. The simplified method provides an efficient and 
easy way to correlate the solubilities of high-molecular-weight solids in supercritical 
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fluids; since no critical properties of the solute are required. 

The unified approach adopted in the rigorous EOS modeling of solid-supercritical 
fluid equilibria integrates different two-parameter cubic EOS with different mixing 
rules in a cohesive manner. It should be especially advantageous in testing the capa- 
bilities of various multiparameter cubic EOS/mixing rule combinations in correlating 
and predicting solid-supercritical fluid equilibria. 
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Nomenclature 


a mixture cohesive energy parameter in cubic EOS 
a,\2 unlike cohesive energy parameter in the mixing rules 
b mixture repulsive energy parameter in cubic EOS 
/> ia unlike repulsive energy parameter in the mixing rules 
ff fugaeity of component i in fluid phase 
ff fugaeity of component i in solid phase 
A: 1 2 binary interaction parameter in the mixing rule for a 
l V 2 binary interaction parameter in the mixing rule for b 
Mi molecular mass of component i 
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I 2 vapor pressure of the solid 
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T temperature 
T e critical temperature 
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v molar volume 
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P a parameter in simplified cubic EOS model = b 2 /by 
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a cubic EOS dependent parameter listed in Table 3.1 
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Chapter 1 


Introduction 


Supercritical fluid extraction technology has been under accelerated development in 
recent years. A myraid of potential application of supercritical fluid extraction (SFE) 
has emerged in disciplines such as biotechnology, food processing , environmental 
control and pharmaceutical and petroleum industry (Brennecke and Eckert, 1989; 
Johnston and Peck, 1989; Larson and King, 1986). Several fundamental aspects of 
the problem have to be looked into and solved before the concept reaches the indus- 
trial arena. These fundamental aspects consist of thermodynamics; including solu- 
bility measurements and modeling; transport phenomena; including measurements 
and modeling of transport properties (viscosity, thermal conductivity and diffusion 
coefficient) and transport coefficients (heat and mass transfer coefficients). 

In many of the applications, especially in biotechnology and food industry, a super- 
critical fluid is used to extract or dissolve high-molecular-weight (HMW) compounds. 
These are usually solid at room conditions and the extent to which they dissolve in the 
supercritical fluid is small. The design of any particular SEE process is based upon 
the knowledge of the solubilities of the low-volatile components in the supercritical 
solvents. One of the major factors limiting the commercial success of SFE is the lack 
of extensive, reliable solubility data for the efficient design of separation units or the 
development of new methods for correlating and predicting supercritical fluid phase 
equilibria. 

Correlation models for the solubilities of solids in supercritical fluids have been 
presented by many investigators, including the application of cubic equations of state 
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(EOS). The key point of the EOS method is the choice of a proper mixing rule to 
determine the mixture parameters. Usually the conventional van der Waals (vdW) 
mixing rules are used with empirically adjusted binary interaction parameters. The 
binary parameters cannot be correlated, and even with the optimal parameters, ap- 
preciable errors still exist for various cubic EOS (Haselow et al.. 1986). This failure 
may be a result of the uncertainty in either the pure compound properties or in the 
solubility data. The experimental solubility data have an error which is often quite 
significant. For example, Johnston et al. (1982) report an experimental error higher 
than 40% for the lower end (below 10~° mole fraction) (if solubility data. In such 
cases, increased correlation errors are expected. 

Recent work in the use of cubic EOS to represent solid-supercritical fluid equilibria 
has focused on the need to develop new mixing rules for highly asymmetric mixtures 
which comprise the majority of the solid-supercritical fluid systems of interest. How- 
ever, the comparison between some recent correlative models using cubic EOS with 
new mixing rules, such as those of Sheng et al. (1992) and Spiliotis et al. (1994), 
with the standard one-fluid vdW mixing rules suggests that there is no significant 
advantage of one method over the other. Any difference is usually within the uncer- 
tainty involved in measuring small solubility values. When dealing with HMW solids, 
therefore, simpler models are usually of more practical interest since the solubility is 
small. 

The lack of reliable critical data for many solid solutes of interest, especially 
complex organic molecules, poses another problem when using cubic EOS to model 
solid-supercritical fluid systems. Since these are needed for conventional EOS model- 
ing, they are commonly estimated using conventional method. Then, the EOS-based 
model has to be fine-tuned by including binary interaction parameters (BIPs) to im- 
prove the quality of the fit. The BIPs so calculated will be highly dependent on 
the method used to estimate the solute’s critical properties. The use of estimation 
methods puts an additional burden on the BIPs since they must correct not only 
the inadequacies of the EOS, but also the uncertainties introduced by the esima- 
tion methods. Recently Estevez et al. (1994) proposed a simplified model using the 
Peng-Robinson (PR) EOS with standard vdW mixing rules to correlate the solu- 
bility of solids in supercritical fluids. In this method the EOS parameters for the 
solute are taken as adjustable parameters rather than attempting to estimate them 
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based on critical properties. The proposed method is considerably simpler than the 
conventional (rigorous) EOS approach and is especially advantageous for correlating 
solubilty of HMW solids in supercritical fluids, since it does not require the values of 
critical temperature or critical pressure of the solute. 

In the present work, several combinations of cubic EOS and mixing rules have 
been studied, computing binary interaction parameters (BIPs) and absolute average 
errors in solubility for binary systems of 14 solid solutes and 5 supercritical solvents at 
different temperatures. Equations of state chosen for this study were Redlich-Kwong 
(RK), Soave (SRK) and Peng- Robinson (PR), whereas the mixing rules chosen were 
quadratic (van dor Waals) and cube -root, using a single BIP, A: 1 2 , or two BIPs, k l2 
and 1 12 ' As expected, the introduction of the second BIP, l l2 , greatly improves the 
accuracy of the calculated solubility values. From the overall percent average errors, 
we can conclude that the cube-root mixing rule performs as well as the conventional 
vdW mixing rule in representing solid-supercritical fluid equilibria and, in general, 
no distinction can be made among the EOS. However, differences among the EOS are 
found when the errors for individual systems are compared. The results of correlation 
using two BIPs are highly satisfactory, considering the uncertainty in experimental 
measurements of low solubility values as well as the uncertainty of the T c , P c and u> 
values of the solute. Promising results have been obtaineed also using the simplified 
model of Estevez et al. (1994) with each of three EOS (RK, SRK and PR) studied. 
Optimum solute-to-solvent parameter ratios a and ft and percent average errors have 
been computed for binary systems of 15 solid solutes and 5 solvents at different 
temperatures. The simplified method provides an efficient and easy way to corelate 
I1MYV solids solubilities in absence of the critical properties of the solute. 

The unified approach adopted in the rigorous EOS modeling of solid-supercritical 
fluid equilibria integrates different two-parameter cubic EOS with different mixing 
rules in a cohesive manner. It should be especially advantageous in testing the capa- 
bilities of various multiparameter cubic EOS/mixing rule combinations in correlating 
and predicting solid-supercritical fluid equilibria. 




Chapter 2 


Literature Survey 


Supercritical fluid extraction (SFE) is a relatively new separation technique that has 
received much attention in recent years (McHugh and Krukonis, 1994). Conventional 
distillation methods involve higher temperatures which result in decomposition of 
heat, -labile substances. Conventional extraction methods involve chemical solvents 
which, even in residual quantities, are of concern from the standpoint of toxicity to 
consumers of the end products. Supercritical extraction, however, is performed at 
relatively low temperatures using solvents which can be easily separated from the 
end products by simple physical means and is a new alternative to distillation or 
liquid extraction (Parkinson and Johnson, 1989). Among the supercritical solvents, 
carbon dioxide is usually preferred, since it is nontoxic, nonflammable, environmen- 
tally acceptable and relatively inexpensive, and its low critical temperature (31.2 °C) 
allows extraction at relatively low temperatures convenient for heat-labile substances. 
SFE exploits the advantage of a solubility enhancement of 10 3 -10 8 in supercritical flu- 
ids and the phenomenon of enhanced solubility in a near-critical fluid of substances 
which are solid at the temperature of the system is well-documented (Robin and 
Vodar, 1953 : Paulaitis et al, 1983: Johnston and Penninger, 1989; Bruno and Ely, 
1991 ). This enhanced solubility is related primarily to the ”liquid-like" density of 
the fluid which promotes strong attractive forces. Additional attractive features of 
SFE are that the low viscosity of the supercritical fluid combined with high solute 
diffusivities result in superior mass transfer characteristics. Furthermore, significant 
density gradients between the particle-fluid interface and the bulk supercritical fluid, 
which are caused by buoyant forces, result in enhanced mass transfer. This enhanced 
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mass transfer is in contrast, to the traditional liquid solvent extraction, which gener- 
ally takes several hours to perform, requires relatively large amounts of solvent, and 
may result in incomplete recovery. In comparison to the conventional processes, SFE 
offers considerable flexibility for an extractive separation through controlling pres- 
sure, temperature, and choice of solvent (Todd and Elgin, 1955). Its main advantage 
over the conventional ones is that the dissolved extract may be completely separated 
from the supercritical fluid simply by decreasing the fluid density, thus its dissolv- 
ing power. This may be accomplished by expansion to subcritical pressure or by a 
relatively small increase in temperature. The energy requirements may thus be re- 
duced significant ly for volatilization and separation of the heavy component (Kohn 
and Savage;, 1979). SFE has already been applied in a number of processes such as 
decaffeination of coffee, production of enzymes and other pharmaceutical substances. 
SFE is indeed very tempting as a separation process in the presence of thermosensitive 
substances. Recent studies have demonstrated the use of supercritical fluids also in 
environmental control, e.g. for extraction and removal of toxic organic contaminants 
from hazardous wastes and various solids such as sediments and fly ash (Brady et ah, 
1987; Paulaitis et al., 1987; Dooley et al., 1987, 1990). Extraction and recovery of 
aromatic hydrocarbons from contaminated soils and waters by means of supercritical 
fluids have been investigated by Hawthorne and Miller (1987), Roop et al. (1989), 
and Akgerman et al. (1992). 

Although the number of applications using supercritical fluids is increasing con- 
stantly, their commercial success is still limited. Solubility dependence of the ex- 
tracted components on pressure and temperature forms the basis of design for condi- 
tions in the extractor and separator. One of the major factors limiting the commercial 
success of SFE is the lack of extensive, reliable solubility data for the ellic.ient design 
of separation units or the development of new methods for correlating and predict- 
ing high-pressure and supercritical fluid phase equilibria. Tin; development of new 
applications for supercritical fluids increasingly implies in-depth knowledge of solute- 
solvent systems, the complexity of which is due to the highly non-ideal behaviour 
of the solutions. The appropriate mathematical modeling of the thermodynamic be- 
haviour of the solute-solvent system is important for effective design of the extraction 
and separation units of a SFE process. Some current treatments attempt empirical 
relations (Chrast.il, 1982; Adaelii and Lu, 1983; Schmitt and Reid, 1985; Kramer and 
Thodos, 1988; Iwai et ah, 1991; Sakaki, 1992; Wen g and Loo, 1992; and Pulitzer et 
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al., 1993) or do modeling that focuses on the virial equation of state (Ewald et al., 
1953; Hinckley and Reid, 1964; Chueh and Prausnitz, 1967; Quiram et al., 1994) and 
cubic equations of state (Dobbs and Johnston, 1987; Dobbs et al., 1987; Kosal ct al., 
1992, Cabellero and Estevez, 1991, 1992; Caballero et al., 1992; Mukhopadhyay and 
Rao, 1993; Yu et al., 1995; Coutsikos et al., 1995; Chen et al., 1995). Even more 
complex equations of state (EOS) have been used to represent fluid-phase nonide- 
ality in correlating solid-fluid equilibria (Finck et al., 1992; Bamberger et al., 1994; 
Bamberger and Maurer, 1994; Deak and Kemeny, 1994; Barna et al., 1994). Since 
experimental measurements, at pressures as high as those involved in SFE, are rather 
difficult, a model for prediction of solubilities is very desirable for design purposes. 
However, recent attempts at the prediction of solubility of aromatic solids in super- 
critical CO 2 with a group-contribution equation of state (Bamberger and Maurer, 
1994; Bamberger et al., 1994) or with the use of predictive excess Gibbs free energy 
(G E ) models in the attractive term parameter of cubic EOS (Spiliotis et al., 1994) - 
the so called predictive EOS/G e models - has not met with much success; the errors 
being typically in the range of 20-50%. 

The empirical correlations, which are often described as linear relations, rarely 
show the behaviour of enhanced solubility of solids in supercritical fluids over wide 
ranges of pressure (Knez and Steiner, 1992). Recently, a linear relationship has 
been developed (Wang and Tavlarides, 1994) based on a dilute two-region solution 
theory. This relationship has the capability to quantitatively describe the solubility 
behaviour of a heavy solute in a compressed gaseous solution. Yu et al. (1995) 
obtained a good correlation of the linear equation to the experimental solubility data 
of polychlorinated biphenyls (PCBs) with the average deviations ranging from 11 to 
19%. A problem often faced with the correlation - and of course prediction - of 
solid solubilities in supercritical gases is the lack of vapor pressures of solutes (e.g. 
solid n-alkanes). In such cases, a solution model, based on the assumption that the 
supercritical phase is dense enough to behave like a (hypothetical) liquid (expanded 
liquid treatment), that describes solid-expanded liquid equilibria has been applied 
to correlate the solubility data (Iwai et al., 1992; Spiliotis et al., 1994). It has the 
advantage that only the molting point of the solute and its enthalpy of melting, which 
are usually available, are needed. In general, however, cubic EOS or even EOS derived 
from statistical mechanics and lattice gas models have been applied in the correlation 
of solubility data with good results. Since EOS derived from lattice gas models are 
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more complicated than a simple cubic EOS (e.g. Johnston et ah. 1982; Balbuena et 
ah, 1991). we focus our attention upon the cubic ones. 

Cubic EOS e.g. Redlich and Kwong (1949), Soave (1972) and Peng and Robinson 
(1976) have been widely employed to correlate solubilities of solids in supercritical 
fluids. These EOS and their temperature-dependent functions were originally derived 
to decribe YLE systems. It has been suggested that these EOS may not be appro- 
priate to describe systems in which solid solutes are involved. As pointed out by 
Schmitt and Reid (1986), the use of these equations to deal with solid solutes implic- 
itly assumes the solid phase to be a subcooled liquid, and as such, the integration step 
used when equating fugacities proceeds essentially from a vapor to a liquid phase. In 
spite of this, these equations have been found to correlate binary solid-supercritical 
fluid experimental data quite well in region away from criticality, especially data for 
non-polar symmetric systems. Lee et al. (1988) examined the Redlich-Kwong (RK), 
the Soave (SRK). and the Peng-Robinson (PR) EOS to evaluate their capability for 
correlating solubility data of solids in compressed gases. They correlated available 
experimental data for 35 binary systems using both one and two binary interaction 
parameters (BIPs) in the conventional van der Waals (vdW) mixing and combin- 
ing rules to express the composition dependence of EOS parameters. Results were 
reported in terms of one average BIP or one pair of average BIPs applicable to a 
range of experimental temperatures and pressures. This way of correlating BIP val- 
ues might not be appropriate if these values are intended for predictive purposes. It 
is well-known that at least the BIP used in the cohesive energy parameter is tem- 
perature dependent and this dependency is strongest for those solid compounds with 
low melting points and/or solubilities above the mole fraction of 0.01 at the highest 
temperatures and pressures (Bartle et al., 1991). Recent work in the use of cubic EOS 
to represent solid-supercritical fluid equilibria has focused on the need to develop new 
mixing rules for highly asymmetric mixtures which comprise the majority of the solid- 
supercritical fluid systems of interest and to cope with the limitations imposed by the 
use of conventional corresponding states theory based on critical properties (Rao and 
Mukhopadhavay, 1988; Soave, 1991: Spiliotis et al., 1994). Rao and Mukhopadhayav 
(1988) proposed a covolume dependent mixing rule for the cohesive (attraction) en- 
ergy parameter of the PR and SRK EOS. Results show a reduction in the sensitivity 
of the calculated solubility values to the BIP. Spiliotis et al. (1994) employed the cu- 
bic t-mPR EOS (Magoulas and Tassios, 1990) to correlate the solubility of aromatic 
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solids in CO 2 using the conventional one-fluid vdW mixing rules with either one or 
two BIPs and also the EOS/G e model based on the UNIQUAC G E expression. In 
addition, an EOS/G E model which coupled the SRK EOS with UNIQUAC was also 
investigated. Very good results were obtained with the t-mPR EOS using either the 
EOS/G E model or two BIP conventional mixing rules, with the latter possessing a 
slight advantage. The comparison between correlative models, such as those of Sheng 
et al. (1992) (modified EOS/G E ), Johnston et al. (1982) (augmented vdW), and the 
t-mPR EOS models suggests that there is no significant advantage of one method over 
the other. Any difference is usually within the uncertainty of the experimental data 
and cannot be considered significant. Chen et al. (1995) proposed a modified mixing 
model with a volume correction term to calculate the solid solubilities in supercritical 
C0 2 and cosolvents using the PR EOS coupled with the modified UNIFAC group 
contribution G E model (Larsen et al., 1987; Sheng et al., 1989). This volume correc- 
tion term is well-correlated as a function of the molar volume of the solid molecule. 
This mixing model can be applied to various EOS with the same type of correlations 
for the volume correction constants (Chen et al., 1994). With this mixing model, 
solid solubilities in supercritical fluids can be predicted with acceptable accuracy. 

The lack of reliable critical data for many solid solutes of interest, especially 
complex organic molecules, poses another problem when using cubic EOS to model 
solid-supercritical fluid systems. The usual approach has been to use several estima- 
tion methods (Reid et al., 1987). The use of estimation methods puts an additional 
burden on the BIPs since they must correct not only the inadequacies of the EOS, 
but also the uncertainties introduced by the estimation methods. To avoid the use of 
estimated properties, Schmitt and Reid (198G) used a modified Peng-Robinson EOS 
in which the pure component parameters for the solid solutes were treated as vari- 
ables, and no BIPs were used in the conventional vdW mixing rules. The modified PR 
EOS with its two adjustable parameters for the solute correlated the solubility data 
better than the classical PR EOS with a single BIP. Recently Estevez et al. (1994) 
proposed a simplified method using the PR EOS with standard vdW mixing rules to 
correlate the solubility of solids in supercritical fluids. It is especially advantageous 
for correlating solubility of high-molecular- weight solids, since it does not require the 
critical properties of solutes. 

Recent interest in supercritical fluid research has centered on the behaviour of 
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supercritical fluid mixtures near critical points. The formation of macroscopic clusters 
of solvent about the solute has been observed to occur at points of criticality resulting 
in an increase in solubility of the solute in the solvent (Kim and Johnston, 1987). 
Conventional cubic EOS are incapable of predicting this kind of behaviour. More 
complex models have been developed taking into account the clustering phenomenon. 
An excellent review on general modeling of supercritical mixtures, including the use 
of cubic EOS. has been made by Johnston et al. (1989). In spite of their limitations, 
cubic EOS exhibit a good simplicity/accuracy balance which makes them particularly 
attractive for engineering purposes. 



Chapter 3 


Thermodynamic Modeling For 
Solid-Fluid Equilibria 


3.1 Equation of State Modeling of Solubilities 


At solid-fluid phase equilibrium, the condition of equal fugacities in both phases (s,g) 




(3-1) 


is fulfilled for each component i. With the usual assumptions: the solubility of the 
supercritical fluid solvent (1) in the solid phase is negligible and the solid phase is pure; 
the solid phase is incompressible; the fugacity coefficient of the pure solid component 
(2) at saturation conditions, (f)% lt , is taken to be unity since the vapor pressure of solid 
is very low, the fugacity /| is given by 


fi = P 2 ex P 


V s 2 (P - Pi) 


RT 


(3.2) 


where is the vapor pressure of the solid at temperature T, is the solid molar 
volume, and P represents the system pressure. The fugacity of the solute (2) in the 
fluid phase is given by 

fi = V2<t>2P (3-3) 


where y 2 is the solubility (mole fraction) of the solid in the supercritical fluid phase, 
and (f >2 is the fugacity coefficient of the solid in the supercritical phase. 
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where y-2 is the solubility (mole fraction) ol the solid in the supercritical fluid phase, 
and (j ) 2 is the fugacity coefficient of the solid in the supercritical phase. 


brom the condition of equal fugacities for the solid in both phases (Equation 3 . 1 ). 
the solubility of the solid in the fluid phase is 

.. . p ,gyK(P - pp/m 


Experimental pure component vapor pressure values (P|), and the solid volume 
v 2 can he taken from the literature. can be calculated from an EOS for a given 
■set ol pressure and temperature of the specified system. 
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The composition dependence of parameters a and 6 for the mixture is represented 
by the conventional vdW mixing rules. For a binary mixture. 

a — ^liy\ + ^y\y20-\2 + ^22V2 ( 3 - 8 ) 


and 


b — bny\ + 2yiy 2 bi2 + 6222/2 


(3.9) 


where iji is the mole fraction of component i , and an and bn are given by Equations 3.6 
and 3.7, respectively. The combining rules for the cross parameters a i2 and 612 are 


a 12 — \Abl&22(l ~ k 1 . 2 ) 


(3.10) 


and 

6 12 = j ^ - lv2 ^ ( 3 - n ) 

where and ly> are the binary interaction parameters. When / i2 is set equal to zero, 
Equation 3.9 simplifies to 


2 

6 = 53 Vi b i = Vi 6 i + 2/262 
2=1 


(3.12) 


A cube-root, mixing rule for the parameter a was recently proposed by Pongsiri 
and Viswauath (1989) and corrected by Caballero and Estevez (1991). It is given by 


‘ i/s - yi ( ;; S / 3 ::^ )1/3 + »» (i - *■>> 


( 3 . 13 ) 


where Mi s the molecular mass of component i. The parameter 6 is calculated using 
the conventional vdW mixing and combining rules. 


3.3 Fugacity Coefficient from Cubic EOS 


The fugacity coefficient fc is given by (Prausnitz et al., 1986): 


In (j)i 


r 00 

' fdP s 

RT 

v 

\dn i/ 

1 

% 

> 

Fc 


dV — In 2 


( 3 . 14 ) 


where Vis the toteJt- volume and z = Pv/RT is the compressibility factor of the mix- 
ture, and rii and rtj are the mole numbers of component i and j, respectively. Equa- 
tions of state which are pressure-explicit can now be used to determine the analytical 
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form of ( dP/dn^TXnj for calculating the fugacity coefficient. However, since EOS 
mixing rules directly relate to mole fractions rather than to mole numbers, deriva- 
tives with respect to mole fractions would greatly facilitate the treatment of various 
mixing rules in cubic EOS for introducing functional modifications in the expres- 
sion for component fugacity coefficients. Szarawara and Gawdzik (1988) developed 
a cohesive technique of generalizing the composition dependence of the component 
fugacity coefficients for the various types of two-parameter cubic EOS which represent 
P = F(T, v. a, b): 


111 Oi 


Pv 

RT~ + RT 


11 Jv Jv° V 


[%-Zyj 


da 

dyj 


f 

Jv 


°° dP 


RT Jv da 
where v° is the molar volume of the ideal gas 


dv + 


' Ob 




db 


ViWi\ r d ? 


RT 


db 


dv (3.15) 


The generalized Schmidt- Wenzel EOS (Equation 3.5) with A > 0. where A = ir — iw 
represents the cubic EOS of Table 3.1 used in the present study. Equation 3.15 may 
be integrated in a general form when Equation 3.5 is used with A > 0. The result is 


In <t>i 


(h_ _I±yjk + ^ - l) - In (z - B) 

-4 [ f(ii-Zyjdj\ fbj -YjjjbA 

BV AL V a J { b J 



(3.16) 


where 


( + \/A 

\ (-U--JX 

); r,-( 2 

V 2 

db . 

* " %’ 


da 

O'i — X j 

oyi 

r — ' _ v — ' da 


(3.17) 

(3.18) 

(3.19) 


(RT) 2 

B = A 

RT 

Further transformations of Equation 3.16 depend on the mixing rules. 
Use of the quadratic vdW mixing rules 

= 


(3.20) 

(3.21) 


a 


(3.22) 
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and 


results in 


^ — £ £ ViVjbij 

* 3 

(3.23) 

“< = ^- = 2 E» a y 

(3.24) 

= z«wr 2a 

(3.25) 

l, = ^ = 

(3.26) 


(3.27) 


In Equations 3.24 and 3.26, a,-j and 6^ are given by Equation 3.10 and 3.11 respec- 
tively. 


If cube-root mixing rule for a is employed (Equation 3.13), we get 


= da = 3a§ 
0 m 


Mj (an M{) 3 Mf (ajjMj ) 3 — — — u ^ n 0 „, \ 

; 772 {v^n a 22 ) i (1 — kij) (1 — 2 Vi) 


(3.28) 


and 


£ w = £ 2 /j 


3a 

dVj 


3a 5 


{Vi ~ Vj) i (afiMi) 3 - Mj 3 {a^MjY 


1\2 
" 3 


+3a3 


El/fcAA 8 

(v^n a jj)~ (1 ” ^ij) { 1 2 ^ + 2/j)} 


(3.29) 


3.4 A Simplified Cubic EOS Model for Solubility 

The simplified model presented below is especially advantageous for correlating sol- 
ubility of high-molecular-weight (HMW) solids in supercritical fluids, since it does 
not require the critical properties of the solute (which are often not available). This 
model has the following key features (Estevez et al., 1994): 
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will be highly dependent oil the method used to estimate the solute’s critical 
properties). 

• Since solubilities are usually low, infinite-dilution compressibility factor is used 
to compute the infinite-dilution fugacity coefficient of the solute in the super- 
critical fluid phase. 


The first feature renders the use of critical property estimation methods needless. The 
second feature tremendously simplifies the computational aspect of the model, since 
it eliminates the need for an iterative scheme. Using the infinite-dilution concept 
renders the compressibility factor and the fugacity coefficient of the solute in the 
supercritical fluid-phase independent of the solubility of the solute. The conventional 
iterative procedure is thus substituted by a once-through calculation. 


The modifications or simplifications herein can be used with any EOS. In this 
work, the RK, SR.K and PR EOS have been used to illustrate the method. For these 
EOS, Equation 3.16 gives 


In fa = j {z - 1) - In (z - B) - ■—= 


2(yiai2 + y-zai) b 2 

In 

' z — r 2 B' 

a b 


Iz-nBl 


(3.30) 


As can be seen, y 2 appears explicitly and implicitly through z which is a function of 
composition in the supercritical phase. For the fluid parameters, a and b, standard 
vdW mixing rules are used (Equations 3.8 and 3.12). 


For the infinite dilution limit (y 2 = 0, y x = 1, a = a x and b — bi), Equation 3.30 
reduces to 


ln<*r = - 1) - K*” - B > - iTK 


2a\2 b 2 
a-i h 


In 


r 2 B 


(3.31) 


Lz°° -nBJ 

Using the standard combination rule for a i2 (Equation 3.30 with k X2 — 0), Equation 3.31 
can be written as 


ln#°“f( 


^2 / „O0 


1) - ln(z°° - B) - 


A 


2 I°1 -h 

bV a L V 


111 


r 2 B] 


r x B\ 


(3.32) 


a = 




(3.33) 

(3.34) 


Dcfinining parameters a and (5 as 
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Definining parameters a and (5 as 


(3.33) 

0 = r (3.34) 

h 

the final equation for the infinite-dilution fugacity coefficient of the solute in the 
supercritical fluid phase is given by 



In <f>? = (5{z°° - 1) - ln(z°° - B) - 


A{2a — /3) 
£v/A 


~z°° — r 2 B' 
_z°° — r\B . 


(3.35) 


In all the equations for Iik/j^ 0 (Equations 3.31, 3.32, 3.35), A and B are defined as 
infinite dilution limits of Equations 3.20 and 3.21, respectively. Reduced properties 
are normalized with respect to the solvent critical properties. 



Chapter 4 


Results and Discussion 


4.1 Physical Properties 


The physical properties ( M,T C ,P C and to) of the supercritical solvents are presented 
in Table 4.1. The solid solutes used in this study are presented in Table 4.2, along 
with their critical properties (T C) P c ), molecular weights, acentric factors and molar 
volumes. The constants of the vapor pressure relation for the solid solutes are also 
presented in Tahiti 4.2. Table 4.3 presents the data for the sublimation pressures of 
PCBs at the temperatures of interest. The references from which the data were taken 
are also shown in the tables. 


Table 4.1: Physical Properties of the Supercritical Solvents 




T c 

Pc 





Solvent 

M.W. 

(K) 

(bar) 

LO 

References 

C0 2 

44.010 

73.8 

304.1 

0.239 

Reid 

et, al., 

1987 

C 2 Ho 

30.070 

48.8 

305.4 

0.099 

Reid 

et al., 

1987 

0 2 ll, 

28.054 

50.4 

282.4 

0.089 

Reid 

et al.. 

1987 

OC1F;, 

104.459 

38.7 

302.0 

0.198 

Reid 

et al. 

1987 

CIIF ;i 

70.013 

48.0 

299.3 

0.200 

Reid 

et ah. 

, 1987 
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Table 4.3: Sublimation Pressures of Polychlorinated Biphenyls 


Solute 

P snt (bar) 

References 

308.1 K 

313.1 K 

323.1 I< 

4, 4'- DOB 

NA 

1 .986 

6.594 

Yu et al., 1995 

2,2 / ,3,3',4,4'-HCB 

NA 

0.2845 

1.0104 

Yu et ah, 1995 

2,3',4 , ,5-TCB 

0.2359 

0.4546 

1.5886 

Yu et ab, 1.995 


4.2 Data Reduction and Parameter Estimation Cri- 
teria 


In the rigorous EOS modeling of the solubilities, the binary alteration parameters 
kjj and Ijj which account for the solute-solvent energetic interactions and size and 
shape dilfercnces, respectively were determined by a regression method minimizing 
an objective function which measures the deviations between experimental ( y 2 ,ex P ) 
and calculated ( 7 / 2 ,c<a) values of the solute solubility. The objective function (OF) 
used is the root-mean-square (mis) of the absolute deviations, which is defined as 


OF = 




(4.1) 


where N is the number of experimental data points. Rosenbroek’s optimization algo- 
rithm (Beveridge k Schechter, 1970) is used to minimize OF. The quality of the fitting 
of the EOS to the experimental solubility data is measured in terms of a percentage 
average error, AE (%), defined as 



,/EL 

,-) 2 

V N 



U2 ,exp 


OF 

'!J2,t'xp 


x 100 = N— 

E 


x 100 
OF 

N 

j=l !J2,c.xp 


x 100 


(4.2) 


The listing of the computer program used to obtain the optimum values of the BIPs 
appears in Appendix A. 


The parameters a and of the simplified model were determined by minimizing 
an objective function which measures the relative deviations between experimental 
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and calculated values of the solubility. The objective function (OF) used is: 


OF = 



2/2, eip 2/2, cai 
2/2, exp 


2 


(4.3) 


The listing of the computer program used to obtain the optimum values of the pa- 
rameters a and (3 appears in Appendix B. 


4.3 Results and Discussion 

Correlation results for binary systems of 14 solutes and 5 solvents at different tem- 
peratures for each of the three EOS tested are given in Tables 4.4 through 4.7. 
Tables 4.4 and 4.5 show the optimum k X2 values and percent average errors when 
the EOS mixture cohesive energy parameter a is expressed by the quadratic (QMR) 
and the cube-root (CMR) mixing rule, respectively; the mixture repulsive energy pa- 
rameter, b, of the EOS being expressed by the linear mixing rule. Our k X2 values for 
the PR-EOS using quadratic mixing rules are in very good agreement with the values 
reported by other researchers. Most recently Yu et al. (1995) reported k X2 values of 
0.0G05 and 0.0582 for the system 4,4'-DCB/C02 at 313.1 K and 323.1 K, respectively. 
Our corresponding values are 0.0604 and 0.0571. Bartle et al. (1991) reported a k X2 
value of 0.090 for the system fluorene/C0 2 at 308.15 I< using a simplified form of 
the PR-EOS for dilute solutions; our calculated value for the same system using the 
conventional PR-EOS is 0.0826. 

Table 4.6 and 4.7 summarize the optimum k 12 and l\ 2 values and percent average 
errors when parameter a of the EOS is expressed by the QMR and CMR, respectively 
and the parameter b is expressed by the quadratic mixing rule . The binary interaction 
parameter k l2 measures the deviation from geometric mean intermolecular attractions 
assumed for the unlike cohesive energy parameter a \ 2 , while the binary interaction 
parameter l l2 measures the deviation from arithmetic mean intermolecular repulsions 
assumed for the unlike repulsive energy parameter b X2 . Both k X2 and l x2 , being purely 
empirical parameters, can be either positive or negative. A physical interpretation has 
been ascribed to negative k x2 values as being indicative of the presence in the mixture 
of specific chemical interactions such as hydrogen bonding (McHugh and Krukonis, 
1982). However, not much significance can be given to such an interpretation for a 
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parameter which is essentially empirical. The values of the BIPs {k\ 2 and l\ 2 ) depend 
oh the accuracy of the physical properties (T c , P c and to) used to evaluate the pure 
component parameters of the EOS, the source of the experimental solubility data, and 
the lorm ol the objective function used to minimize the deviations between calculated 
and experimental values of solubility. No definite trends can be observed in the 
behaviour ol the k\ 2 and l\ 2 values with increasing temperature. For any particular 
system, however, the behaviour of the ^2 values mirrors that of the k i2 values in 
general, that is, if k\ 2 increases with temperature, so does / 12 ; if fc 12 decreases, l\ 2 also 
decreases. 


The results of correlation shown in Table 4.6 and 4.7 can be considered satisfactory, 
considering the uncertainty in experimental measurements oflow solubility values, e.g. 
10-20 % below 10~' 1 mole fraction and increasing to 40 % for the lower end of the 
solubility data (below 1()~ 5 mole fraction) as reported by Johnston el, al. (1982), as 
well as the uncertainty of the T c , P c and to values of solute. Comparison of the overall 
percent average errors of Table 4.4 and 4.5 with those of Table 4.6 and 4.7 shows that, 
as expected, the introducion of the second binary interaction parameter, l l2 , greatly 
improves the accuracy of the calculated solubility values. The calculated solubilities 
of anthracene in supercritical carbon dioxide at 323.15 K using the PR EOS/QMR 
formulation are shown in Figure 4.1. The plots clearly demonstrate that using the two 
BIPs (k\ 2 and l\ 2 ) results in a marked improvement in the fit of calculated solubilities 
with the experimental data in comparison to that resulting from the use of k v2 alone 
or setting k l2 = 0. 

From the overall percent average errors of Tables 4.4 to 4.7, we can conclude 
that the cube-root mixing rule performs as well as the conventional vdW mixing 
rule in representing solid- supercrtical fluid equilibria and, in general no distinction 
can be made among the EOS. However, on closer inspection, differences among the 
EOS are found when the errors for individual systems are compared. Some typical 
results are shown in Figures 4.2- 4.5. The calculated solubilities of benzoic acid in 
fluoroform at 328.25 K using different EOS and a single BIP (k 12 ) in the cube- 
root, mixing rule are shown in Figure 4.2. Figure 4.3 presents the results for the 
system napthalene/fluoroform at 328.15 K using different EOS/CMR formulation 
with optimum k yi and l n values . Figures 4.4 and 4.5 present the results for the 
systems tluorene/ethylene at 343.15 K and benzoic acid/C0 2 at 338 I<, respectively 
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using different EOS with quadratic mixing rule; Figure 4.4 is based on the use of 
optimum ky 2 alone whereas optimum k\ 2 and have been used in generating the 
plots in Figure 4.5. For the systems represented in Figures 4.2 and 4.3, the data are 
best correlated by the SRK EOS; on the other hand, for the systems represented in 
Figures 4.4 and 4.5, the data appear to be better correlated by the PR EOS than 
by the RK or SRK EOS. This observation could indicate that, besides the use of 
BIPs, the functional form of the attractive term of the EOS plays an important role 
in the ability of the equation to account for the effect of nondispersive-type forces 
in solid-supercritical fluid systems. In systems where dispersive forces predominate, 
no significant difference is likely to be found in the ability of the equations to model 
solid-supercritical fluid equilibia. When specific chemical forces are present, however, 
the nature and strength of these forces, as determined by the type of solvent, will 
influence the ability of these equations to model the system. This ability will depend 
on how well the functional form of the attractive term of the EOS can account for 
the effect of these chemical forces. 


For some systems the nature of the supercritical solvent is observed to have a 
marked influence on correlation results (Tables 4.6 and 4.7). For example, in ben- 
zoic acid /solvent systems a considerable increase in percent average error is noticed 
when the supercritical solvent is changed from CO2 to C2H4. A possible explana- 
tion for this observation could be related to the differences in the polarizabilities of 
the two solvents and the effect of this difference on the solute-solvent interaction 
constant a l2 . The polarizability of C2H4 is 42.52 x IQ -25 cm 3 , while that of C0 2 
is 29.11 x 10~ 2,r> cm 3 (Lido, 1990); the ratio of the polarizability of C 2 IU to that of 
C0 2 being approximately 1.5. For the same solute, and upon changing the solvent 
from CO2 to C2II4 at the same temperature, it is therefore to be expected that the 
unlike interaction constant a.12 should be considerably larger for C2H4 than for CO2. 
However a much lower ratio of a l2 in C 2 H 4 to that in C0 2 is predicted using critical 
properties. The increase in percent average error may point at the inadequacy of 
using critical properties to predict successfully the interaction constant in mixtures 
that are highly asymmetric. 
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Figure 4.1: Calculated Solubilities of Anthracene in Supercritical CO- 2 at 323.15 K Using 
the PR EOS/QMR formulation: Role of Binary Interaction Parameters 
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Figure 4.2: Calculated Solubilities of Benzoic Acid in Fluoroform at 328.25 K Using 
Different EOS/CMR Formulation with Optimum k v2 



0 . 000-1 1 1 1 — 1 — 1 — 1 — *■ — - 1 — 1 

0 50 100 150 200 250 300 350 400 450 500 


Pressure, [bar] 


Figure 4.3: Calculated Solubilities of Naphthalene in Fluoroform at 328.25 K Using Dif- 
ferent EOS/CMR Formulation with Optimum fc 12 & hi 
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Figure 4.4: Calculated Solubilities of Fluorene in Ethylene at 343.15 K Using Different 
EOS/QMR Formulation with Optimum k\ 2 
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Figure 4.5: Calculated Solubilities of Benzoic Acid in C0 2 at 338 K Using Different 
EOS/QMR Formulation with Optimum k x2 &. hi 
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Figure 1.6: Calculated Solubilities of Benzoic Acid in Supercritical COo at Different 
Temperatures Using the Simplified Model with the PR EOS 

Due to nonavailability of the critical properties for the solutes acridine and 1,4- 
naphthoquinone, the rigrous EOS modeling was not done for systems made up of 
these solutes. Table 4.8 presents optimum solute-to-solvent parameter ratios a and 
ij and percent average errors for binary systems of 15 solutes (including acridine and 
1,4-naphthoquinone) and 5 solvents at different temperature for each of three EOS 
using the simplified method. In general, both parameters decrease with increasing 
temperatures . The rate of change of ft with temperature is marginally higher (in 
absolute value) than the corresponding value for a. Comparison of overall percent 
average errors shows that no distinction can be made among the EOS. To illustrate 
typical prediction capabilities of the simplified model, four systems have been ran- 
domly chosen, and their solubility curves at different temperatures calculated using 
the PR EOS are shown in Figures 4.6 through 4.9. Also included in these Figures 
are the experimental data. Qualitative observation of the curves indicate that in all 
cases the simplified model gives fairly accurate predictions of solubility. 
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Figure 4.7: Calculated Solubilities of Phenanthrene in Supercritical CO? at Different 
Temperatures Using the Simplified Model with the PR EOS 
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Figure 4.8: Calculated Solubilities of Acridine in Supercritical C0 2 at Different Temper- 
atures Using the Simplified Model with the PR EOS 
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Figure 4.9: Calculated Solubilities of 1,4-Naphthaquinone in Supercritical Ethane at Dif- 
ferent Temperatures Using the Simplified Model with the PR EOS 





Table 4.4: EOS Modeling of Solubility: Optimum k \2 Using the 
Quadratic Mixing Rule (QMR) 
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Table 4.5: EOS Modeling of Solubility: Optimum k y > Using the 
Cube Root Mixing Rule (CMR) 



Benzoic Acid 001F ;J 418. 15 (i 0.1570 1.14 0.5282 15.78 0.1707 8.8(i Schmitt and Reid (1980) 

Benzoic Acid CC1F 3 328.15 6 0.4709 8.45 0.5233 11.84 0.4812 3.80 Schmitt and Reid (1986) 

Benzoic Acid CIIF.s 318.25 5 0.3749 8.43 0.4494 25.12 0.4001 16.18 Schmitt and Reid (1986) 

Benzoic Acid 011F 3 428.25 5 0.4900 1.21 0,1401 21.01 0,1047 14,10 Schmitt and Reid (1980) 
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Naphthalene C»ll (i 308.15 6 0,1501 12.62 0.5-171 22.18 0.5120 18.77 .lulm.sl.oii el al. (1982) 

Naphthalene C 2 II fl 318.15 8 0.4783 10.38 0.5753 23.26 0.5397 19.48 Johnston et al. (1982) 

Naphthale.no 0,1 1 (i 328.15 13 0.5329 12.89 0.6813 33.12 0.6399 28.71 .lulmstun ol. al. ( 1982) 

Naphthalene C2H4 298.15 10 0.3623 15.10 0.4500 29.58 0.4131 26.17 Johnston and Eckert (1981) 

continued on next page 
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Table 4.6: EOS Modeling of Solubility: Optimum k V i & l 
Using the Quadratic Mixing Rule (QMR) 
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Naphthalene C 2 l! (i 1108. If) (i -0.3813 -0.1:182 (i.01 -0.1)0(17 -0.1927 0.07 -0.1107 -0.1531 0.59 

Naphthalene C 2 H G 318.15 8 -0.3380 -0.0806 3.38 -0.0739 -0.2127 5.43 -0.1194 -0.1495 4.76 

Naphthalene C 2 H 6 328.15 13 -0.2873 0.0068 2.55 -0.0295 -0.0989 5.00 -0.0010 -0.0282 5.02 

Naphthalene C 2 U.i 298.15 10 -0.3017 -0.0051 2.91 -0.1200 -0.2980 0.01 -0.1628 -0.2238 5.47 

continued on next page 
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Overall Average Error (%) 7-li.i 7M2 7.28 



Table 4.7: EOS Modeling of Solubility: Optimum k 
Using the Cub t Root Mixing Rule (CMR) 
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continued Grom previous page 

Equation of State j S-R-K j P-R j R-K 

Solute Solvent j Temp. N j ky ly j AE(%) j ky j ly | AE(%) j ky ly j AE(%) 

Phenanthrene C 2 H 6 313.00 4 j 0.1634 -0.0155 . 3.44 , 0.2214 , -0.1756 103 0.2003 -0.0876 121 

Phenanthrene C 2 H 6 333.00 7 i 0.2031 -0.0283 5.11 0.2190 -0.2277 8.21 0.2036 I -0.1286 7.95 
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Table 4.8: EOS Modeling of Solubility: Optimum Solute to 
Solvent Parameter Ratios a & j3 in the Simplified Method 
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4.4 Conclusions 

Several combinations of cubic equations of state and mixing rules have been stud- 
ied and optimum binary interaction parameters are computed for modeling solid- 
supercritical fluid equilibria. A marked improvement is observed in the accuracy of 
the calculated solubilities when the second BIP, l 12 , is included in the analysis. How- 
ever the quadratic and cube-root mixing rules show no substantial difference when 
the overall percent average errors are compared. Furthermore, the three EOS tested 
also give similar results in terms of overall percent average errors. However, differ- 
ences among the EOS are found when the errors for individual systems are compared. 
The following recommendation emerges from these conclusions: any of the three EOS 
along with any of the two mixing rules tested can be used, as long as two binary 
interaction parameters are included, one for the cohesive energy parameter, a, and 
one for the repulsive energy parameter, b. From a practical point of view, the simplest 
EOS and mixing rule may be preferred. 

Promising results have been obtained also with a model which is considerably 
simpler than the rigrous EOS approach since no iteration is required to find the solute 
fugacity coefficient. Each of the three EOS studied gave fairly accurate predictions 
of solubility. The simplified model should prove to be an efficient and easy way to 
correlate high-molecular-weight solids solubilities, since it does not require the critical 
properties of the solute. 

Lastly, the unified approach adopted in the present work to integrate different two- 
parameter cubic EOS with different mixing rules should be especially advantageous in 
the case of multiparameter cubic equations of state (three and four parameters) while 
deriving the composition dependence of the fugacity coefficient of a component in a 
mixture. There are many established cubic equations of state (Mika, 1989; Lawal and 
van der Laan, 1994) whose utility has not been extensively tested because of the lack 
of a cohesive technique of generalizing the composition dependence of the component 
fugacities for different combinations of cubic EOS and mixing rules. It is felt that 
such a unified approach should be of great utility in testing the capabilities of various 
EOS/mixing rule combinations in correlating and predicting solid-supercritical fluid 
equilibria. 
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Appendix A 


Program Listing of the Rigorous 
Model 


#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include "rosen.h" 

#def ine R 83.14 
#define ITER 100 
#def ine MAX_ITER 20 
#def ine TOL 1.0e-5 
#def ine TOLRNS 1.0e-10 
#define N 20 
#define Size 100 
#define sq(x) ((x)*(x)) 


/* . . . Input Variables . . . */ 

char *fnl=" . . /data/ethyln" , *fn2=" . . /data/2_3dmn" , *dat , 

*Pydat = " , . /data/c2h4_3d2" ; 

int opt,mr,eos; 

int i,u,w,N_exp; 

double mi,tcl,pcl,accentricl; 

double mi , t c2 , pc2 , accentr ic2 , pA2 , pB2 , pC2 , Vs , p_vap ; 
double t , P [N] , y2_exp [N] ; 
double omega^a, omega^b ; 
double CONVRG =0.0 ; 


char *fn[3»{ 
"vdW" » 
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"SRK" , 

"PR", 

"RK">; 

FILE *fp,*fq,*fr,*fd,*sol_fp; 


/* • • . Calculated Variables . . . */ 
double al ,a2,bl ,b2,B.l ,B_2; 


/* solmain.c */ 


#include <stdio.h> 

#include <math.h> 

#include "rosen.h" 

#include "solblty . c" 

#include "rosen.c" 

#include "sol .input . c" 

#include "eql.c" 

#include "prdeql.c" 

extern double avg.error (double) ; 
char *sol_fn="sol . out" ; 

void main() 

{ 

int i ,k=0 ,mr_max=2 ; 
double dif , pred.value ; 
extern int no .variables ; 
int variable; 
extern char *dat; 

extern double Orosenbrock) () , temp.vector [] [MAX.VARIABLE] ; 

double Avg.Err [32] ; 

double (*fp.prdeql) () ; 

char *mr_name ; 

char *prsnt_mr []={ 

"QMR.&.QMR” , 

"CMR.&.QMR" , 

"QMRJcJLMR", 

"CMR.&.LMR" , 

"Wong.&.Sand" }; 

char *table_fn[] ={ 

" lp.QMR" , 

" lp.CMR" , 

"2p.QMR H , 

"2p.CMR" , 

"lp.QLMR" , 

" lp.CLMR” >; 


FILE *reslJCp, *res2„fp, *table.fp[7] ; 
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it ( (sol_fp=f open(sol_fn, "at 1 ') ) == NULL ){ 
fprintf ( st derr, ’’Error opening file: , /.s\n , ‘, 
sol„fn) ; 
exit (1 ) ; 

} 


f or (i=0 ; i<=3 ; i++) { 

if ( (table_fp [i] =f open(table_fn[i] , u at”) ) = 
NULL ){ 

fprintf ( st derr,"Error opening file: 7.s\n H , 
table_fp[i]) ; 
exit (1) ; 

> 

> 


tpy „ input () ; 
system^ input () ; 


fprintf (sol_fp/ f \n\t*** SYSTEM ***\n"); 
fprintf (sol.f p , " \n*7,s At V/ #s At*y,s*\n n , f nl , f n2 , Pydat) ; 
fprintf (sol_fp , "\nTemperature : */,lf\tNo. of Data points : '/.dXn 1 ' , 
t ,N_exp) ; 

re : 

printf ("\nEnter the no of Variables 
scanf (”*/,d n ,&no_variables) ; 

if (no_variables){ 

/*if (no ..variables == 2) 
mrjnax = 2 ;*/ 

fprintf (sol__fp , "\n\nThe no. of OPTIMIZATION parameters : , 

no_variables) ; 
f or (mr=l ;mr<=mr jmax;mr++) { 
mr_name = prsnt_mr [ — mr] ; 
mr++; 

if (no.variables == 1) { 
fprintf (table_fp[mr-l] , 

n \n'/,s\ty,s\t'/. . 21f \t'/,d\t" , 
fn2,fnl ,t ,N_exp) ; 

} 

else { 

fprintf (table_fp[mrtl] , 

" \n'/, s \ t '/, s \ t % . 21f \ t */*d\ t ' ' , 
fn2,fnl,t ,N_exp) ; 

} 

fprintf (sol.fp, M \n\nMixing Rule for \"a\" & \"b\ M is '/.s" , 
mr jiame) ; 

f or (eos=3 ; eos<=3 ; eos++) { 

Select JEoS(eos) ; 

Calc„Constants(eos) ; 
rosenO ; 

predLvalue 3 avg.,„error(((* 
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rosenbrock) (temp^vector) ) ) ; 

pred.value = predict_eql(mr , 
eos ,temp_ vector) ; 

fprintf (sol^fp, "\n7.s " , 

dat) ; 

if (no^variables == 1) { 
fprintf (sol_fp , 

"\t kij = 7, . 41f \t Avg_Err = 7..21f ", 
temp_vector [0] [0] , . 
pred_ value) ; 
fprintf (table_f p [mr-1] , 

"7. .41f \t7. .21f \t" , 
temp_vector [0] [0] , 
pred_value) ; 


else { 

fprintf (sol_fp , 

"\tkij = 7..41f\tlij = 7. . 41f \t Avg_Err = 7..21f", 
temp_vector [0] [0] , 
temp__vector [0] [1] , 
pred_value) ; 
fprintf (table_fp[mr+l] , 

" 7, . 41f \t 7. . 41f \t 7. . 21f \t " , 
temp„vector [0] [0] , 
temp_vector [0] [1] , 
pred^value) ; 

} 

/* eql(mr ,eos , temp_vector) ; */ 

printf ("\n**AVG. ERROR : 7.1f \n" , 
pred_value) ; 

> 

if (no ..variables == 1) { 
f close (table_.fp[mr-l] ) ; 

} 

else { 

f close (table_f p [mr+1] ) ; 

} 

} 

goto re; 

> 

> 


/* solblty.c */ 


#include "solblty.h 1 
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#include "mix^rule. c" 
#include "fug_SW_EoS . c" 
#include "fugJ?R_WS. c" 


/* Calculation of a,b */ 

void calc_const (eos , tc ,pc ,accentric, t , emega^a, omega_ib,a, 
b) 

int eos; 

double t c , pc , accentric , t , omega_a , omega_b , *a , *b ; 

{ 

double tr; 
double alpha, beta; 
tr = t/tc; 
switch(eos) 

{ 

case 1 : 

*a = omega_a * R*R * tc*tc / pc ; 

*b = omega__b * R * tc / pc ; 
break; 


case 2: 

*a = omega_a * R*R * pow(tc,2.0) / pc ; 
beta = 0.480 + 1 . 574*accentric - 0.176* 
accentric*accentric ; 
alpha = 1.0 + beta * (1.0 - sqrt(tr)) ; 
alpha *= alpha; 

*b = omega_b * R * tc / pc ; 
break; 

case 3: 

*a = omega_a * R*R * tc*tc / pc ; 
beta = 0.37464 + 1 . 54226*accentric - 0.26992* 
accentric*accentric ; 
alpha = 1.0 + beta * (1.0 - sqrt(tr)) ; 
alpha *= alpha; 

*a *- alpha; 

*b = omega_b * R * tc / pc ; 
break; 

case 4; 

*a = omega^a * R*R * pow(tc,2.5) / (pc* 
sqrt(t)) ; 

*b = omega_b * R * tc / pc ; 

break; 

default : 

break; 

} 

> 


/* Schmtt-Wenzel EOS in > z J form */ 
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double f(A,B,z,u,w) 
int u,w; 
double A,B,z; 

return z*z*z + (u*B-B-l . 0) *z*z + (w*B*B - u*B*B 
- u*B + A)*z - (A*B + w*B*B +w*B*B*B) ; 

> 

/* Differential of S-W-Eos */ 

double df (A,B,z,u,w) 
int u,w; 
double A,B,z; 

{ 

return 3.0*z*z + 2.0*(u*B - B - 1.0) *z + (w*B* 

B - u*B*B - u*B + A) ; 

} 

/* Calculation "vap-z" k vap-vol. */ 

double calc_v(u,w,p,t,a,b, guess) 
int u,w; 

double p,t , a, b, guess; 

{ 

double z,diff ,dfz ,fz, v; 
double A,B; 
int i=0 ; 

B = b*p / (R*t) ; 

A = a*p / (R*R*t*t) ; 
z = guess; 
do 
{ 

dfz = df (A,B,z,u,w) ; 

/★if (df z<TQL) 

{ 

fprintf (stderr , "Denom < TOL in newton raphson for Z" ) ; 
exit (3) ; 

}*/ 

fz = f (A,B,z,u,w) ; 
guess = z - fz/dfz; 
diff = fabs(z -guess); 
z = guess; 

}while ( i++<ITER kk diff>T0L) ; 
if (diff>T0L) /*not converged*/ 
return -100; 

/★printf ( ,, \nz« , /.lf " ,z) ;*/ 
v = R*t*z/p; 
return v; 

> 


double solblty(t ,p,fug2) 
double t »p»fug2 ; 
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double pf; 

pf= exp(Vs*(p-p_vap)/ (R*t)) ; 

return p__vap*pf / (p*fug2) ; 

> 

double func(x) 

double x[] [MAX_ VARIABLE] ; 

{ 

int i , j , stat ; 
double y2_guess, dif; 

double guess=l . 0 , y2_cal [N] , dif f [N] , sum=0 . 0 ; 
double v „by„SW , am , bm , f ug2_SW , OF , k , 1 , f ug2_PR ; 
double Q ,D,dJD,d_Q,gamma_i,Ae ; 

B_1 = bl - al/(R*t) ; 

B„2 = b2 - a2/(R*t) ; 
k = x[0] [0] ; 
if (no_variables == 1) 

1 = 0.0 ; 
else 

1 = x[0][l] ; 

for (i=0 ; i<N_exp ; i++) { 

j =0 ; 

y2_cal [i] = y2_exp[i] ; 
do { 

y2__guess= y2_cal[i]; 
if (mr == 5) { 
eval_Ae(y2_guess,P[i] ,& 

Ae ,&gamma__i) ; 
eval_D_Q_wsmr (y2_guess , 

k , Ae , gamma_i , &Q , &D , &d_Q , 

&d_D) ; 

} 

mr_f or_a _b (y2_guess , k , 1 , Q ,D , &am , 

&bm) ; 

v_by«.SW = calc_v(u,w ,P [i] ,t ,am, 
bm, guess) ; 

/★if (!j && v_by_SW <= bm ) { 

printf ("\n bm is greater than v !!!!!") 

fug2„SW = eval_fugJ5W(k,l,u,w,P[i] , 

v„by-SW , t , am , bm, D , Q , d _D , d_Q , y2_guess) ; 
fug2_PR = eval_fug_PR„WS(y2_guess, 

P[i] ,v„by.SW,t,am,bm,Ae,gamma^i, 
k) ; 

y2_cal[i] = solblty(t ,P [i] ,fug2_SW) ; 
dif = f abs (y2_guess-y2_cal [i] ) ; 

} while (j++ < MAX_ITER kk dif > TOL) ; 

> 

for (i*0;i<NL.exp;i++) 

{ 


; break;} ★/ 
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diffCi] = Cy2_cal[i] - y2_exp[i])*(y2_cal[i] 
- y2_exp [i] ) ; 
sum += dif f [i] ; 

} 

OF = sqrt(sum / N_exp) ; 
return OF ; 

} 


double avg_error (OF) 
double OF; 

{ 

int i; 

double sum=0.0; 

for (i=0 ; i<N_exp; i++) 

{ 

sum += y2_exp [i] ; 

> 

return OF*N_exp*100 . 0/sum ; 

} 


void Select _EoS(eos) 

{ 

switch(eos) 

{ 

case 1 : 
dat=f n [0] ; 
break; 
case 2: 
dat=fn[l] ; 
break; 
case 3: 
dat=fn[2] ; 
break; 
case 4: 
dat=fn[3] ; 
break; 
default: 
break; 

} 

if ( (fp=f open(dat , M rt") ) == NULL ) 

{ 

fprintf (stderr , "fopen failure : 7,s\n",dat); 
exit (1) ; 

> 

f scanf (fp , M %d%d%lf # /.lf " ,&u,&w,&:omega„a > &omega_b) ; 
f close(fp) ; 

> 
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void Calc„Constants(eos) 

calc_const (eos , tel ,pcl , accentricl , t , omega_a, omega_b , 
&al,&bl); 

calc_const (eos , t c2 ,pc2 , accentric2 , t , omega_a, omega.b , 
&a2,&b2); 

> 


/* fug_SW_EoS.c */ 


/* Determinent */ 

double eval_D(u,w) 
int u, w; 

{ 

return ( u*u-4*w ) ; 

} 

/* Roots */ 

void eval_rl_r2(u, w,bm,rl ,r2) 
int u,w; 

double bm,*rl,*r2; 

{ 

double D; 

D = eval_D(u,w); 
if (D<0 . ) 

{ 

printf ("Delta is negative!”); 
exit (1) ; 

} 

*rl = ( (-u+sqrt (D) ) /2) *bm; 

*r2 = ((-u-sqrt (D))/2)*bm; 

} 

/* First term in fugacity expression */ 

double term_01(p,v,t) 
double p,v,t ; 

{ 

return ( p*v/(R*t) -1.0 ) ; 

} 

/* Second term in fugacity expresion */ 

double termj)2(t,p) 
double t»p; 

{ 
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double vd; 
vd = R*t/p; 
return ( log(vd) ) ; 

} 

/* Third term in fugacity expression */ 

double term_03(u,w,v,bm,t,am) 
int u,w; 

double v,bm,t,am; 

{ 

double temp,D,rl,r2; 

D=eval_D(u,w) ; 
eval_rl_r2(u,w,bm,&rl,&r2) ; 
if (D==0.0) 

return ( -log(v-bm) - am/ (R*t* (v-rl) ) ) ; 

> 

if (D>0.0) 

{ 

temp = am * log( (v-r2) / (v-rl) ) / (R*t*(rl- 
r2) ) ; 

return (-log(v-bm) - temp) ; 

> 

> 

/* Fourth term in fugacity expression */ 

double term_04(u,w, v,bm) 
double v,bm; 
int u,w; 

{ 

double D,rl,r2; 

D= eval_D(u,w) ; 
eval_rl_r2(u,w,bm,&rl,&r2) ; 
if (D>0 . ) 

{ 

return ( -log((v-r2)/(v-rl)) / (rl-r2) ); 

> 

if (D==0 . ) 

{ 

return (-1) / (v-rl) ; 

> 

> 

/* Fifth term in fugacity expression */ 

double -erm_05(u,w,a,b,v) 
double b,a,v; 
int u,w; 

{ 

double D,rl,r2; 

double tl,t2,t2a,t3a,t3,t4; 

D=eval_D(u,w) • 
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eval_rl_r2(u,w,b,fcrl,&r2) ; 
if (D==0 . ) 

■c 

tl= (R*t)/(v-b); 

t2= a*u/ (2*pow(v-rl,2.0)) ; 

t3= (2*w*a*b + rl*a*u)/(3*pow(v-rl,3.0)) ; 

return tl+t2+t3 ; 

> 

if (D>0 . ) 

{ 

tl= (R*t)/(v-b); 

t2a= v*v + u*v*b + w*b*b; 

t2= a*u*0.5/t2a; 

t3a= (2*w*a*b-a*u*u*b*0.5 )/((rl-r2)*(rl- 
r2)) ; 

t3= t3a*(2.0*v - (rl+r2))/t2a; 

t4= t3a*2.0*log((v-r2)/(v-rl))/(rl-r2) ; 

return tl+t2+t3-t4; 

> 

> 


/* Fifth (a) in fugacity expression when LMR for "b" */ 

double term_05_a_lmr (bi) 
double bi; 

return (bi) ; 

> 

/* Fifth (b) in fugacity expression when LMR for "V */ 

double term_05_b_lmr(bm) • 
double bin; 

{ 

return bm ; 

> 


/* Fifth (a) in fugacity expression when QMR for "b" */ 

double term_05_a_qmr(bi,bm,lij) 
double bi,bm,lij; 

■c 

return (bi+bm)*(1.0 - lij); 

> 

/* Fifth (b) in fugacity expression when QMR for "b’ */ 

double term_05_b_qmr(bm) 
double bm; 

{ 

return ( 2.0*bm ); 
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> 

/* Fourth (b) in fugacity expression when QMR for "a" */ 

double, term_04_b_qmr (am) 
double am; 

■C 

return ( 2.0*am ); 

> 

/* Fourth (a) in fugacity expression when QMR for "a" */ 

double term_04_a_qmr(ai,am,kij) 
double ai , am , ki j ; 

{ 

return 2.0*sqrt (ai*am) *(1 .0 - kij); 

> 

/* Fourth (a) in fugacity expression when CMR for "a" */ 

double term_04_a_cmr (kij ,x2 ,am) 
double kij,x2,am; 

{ 

double xl,tl,t2,t3,t4; 
xl = 1.0 - x2; 

tl= pow(mi, . 3333) *pow(a2*mj , .3333) ; 

t2= pow(mj , . 3333) *pow(al*mi , .3333) ; 

t3= sq(xl*pow(mi , . 3333) + x2*pow(mj , . 3333) ) ; 

t4= pow(sqrt (al*a2) , .3333) * Cl— kij ) * (1 - 2*x2) ; 

return( 3*pow(am, . 6667) * ( C(tl-t2)/t3) +t4) ); 

} 

/* Fourth (b) in fugacity expression when CMR for "a" 

double term_04_b - cmr (kij ,x2,am) 
double kij,x2,am; 

1 

double xl,tl,t2,t3,t4; 
xl = 1.0 - x2; 

tl = pow(mi, .3333) *pow(a2*mj , .3333) ; 
t2=pow(mj , . 3333) *pow (al*mi , .3333) ; 
t3=sq(xl*pow(mi , . 3333) +x2*pow(mj , .3333)) ; 
t4=pow (sqrt (al*a2) , . 3333) * ( 1-ki j ) * ( 1-2* (x2*x2+xl* 
xl)) ; 

return ( 3*pow(am, . 6667) * ( ( (xl-x2)*(tl-t2)/t3) + 
t4) ) ; 

> 

/* Fourth (a)&(b) in fugacity expression when WS-MR */ 
void t erm_04_a„b jwsmr (bm , D , d^D , t5a , t5b , t4a , t4b) 
double bm,D > d ww D»t5a,tSb,*t4a,*t4b ; 

{ 


*t4b 


R*t*(bm*D + D* t5b ) ; 
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*t4a = R * t * ( bm * d_D + D * t5a) ; 
return ; 

> 

/* 

void term_04_a_b_wsmr (D , Q , dJD , d_Q , t4a , t4b) 
double D,Q,d_D,d_Q,*t4a,*t4b : 

double tempi, temp2 ; 

tempi = Q/sq(1.0 - D) ; temp2 = D/(1.0 - D) ; 

*t4b = R * t *( D * tempi - 2 * Q * temp2 ) ; 

*t4a = R * t * ( tempi * d 3 + temp2 * d_Q) ; 

} 

*/ 

/* Fifth (a)&(b) in fugacity expression when WS-MR */ 

void term_05_a„b_wsmr (D, Q ,d_D ? d_Q ,t5a,t5b) 
double D,Q,d_D,d„Q,*t5a,*t5b ; 

{ 

double tempi, temp2 ; 

tempi = Q/sq(1.0 - D) ; 
temp2 = D/(1.0 - D) ; 

*t5b = D * tempi + 2 * temp2 * (Q/D) ; 

*t5a = tempi * d_D + d_Q / (1.0 - D) ; 

} 

void term_04_a_b(mr ,ai,am,x2,kij ,t4a,t4b) 
int mr; 

double ai , am,x2 ,kij , *t4a, *t4b ; 

{ ' 

switch(mr) 

{ 

case 1 : 

*t4a = term„04_a_qmr (ai,am,kij ) ; 

*t4b = term„04_b_qmr (am) ; 

break; 

case 3 : 

*t4a = term_04_a_qmr(ai,am,kij) ; 

*t4b = term_04^b^qmr (am) ; 

break; 

case 2 ; 

*t4a = term^.04^a„cmr(kij ,x2,am) ; 

*t4b * term !M .04™.b..cmr(kij ,x2,am) ; 

break; 

case 4 ; 

*t4a * term_04_a_cmr(kij ,x2,am) ; 

*t4b 58 term.04 a .b w cmr(kij ,x2,am) ; 
break ; 
default : 
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break ; 

> 

} 


void term_05_a_b(mr ,bi ,bm,lij , t5a, t5b) 
int mr; 

double bi,bm,lij *t5b; 

{ 

switch (mr) 

{ 

case 1 : 

*t5a - term_Q5_a_qmr(bi,bm,lij) ; 

*t5b = term_05_b_qmr (bm) ; 

break; 

case 3 : 

*t5a = term_05_a__lmr (bi) ; 

*t5b = term„05„b_lmr(bm) ; 
break ; 
case 2 : 

*t5a = term_05_a_qmr (bi ,bm, lij ) ; 

*t5b = term„05_b_qmr (bm) ; 

break; 

case 4 : 

*t5a = term_05_a_lmr (bi) ; 

*t5b = term_05_b_lmr (bm) ; 

break; 

default : 

break; 

} 

} 

/ * Evaluate f ugacity of a component */ 

double eval_fug_SW(k,l,u,w,p, v,t ,am,bm,D,Q ,d_D ,d_Q,x2) 
int u,w; 

double p,v,t,am,bm,k,l,x2,D,Q, d_D , d_Q ; 

{ 

double tl ,t2 ,t3 ,t4,t5 ,t4a,t4b,t5a,t5b; 

tl = term_01(p, v,t) ; 

t2 = term_02(t,p) ; 

t3 = term_03(u,w,v ,bm,t ,am) ; 

if (mr-=5) { 

term_0S_a Jd jwsmr (D , Q , d_D , d^Q , &t 5a , &t5b) ; 
t erm w .04_a.b„wsmr (bm , D , d_D ,t5a , t5b , &t4a , 

&t4b) ; 

/* 

t erm.04._a Jj^wsmr (D , Q , d^D » d_Q , &t4a , &t4b) ; 

*/ 

} 

else { 

term*04 w ,a„b (mr ,a2,ttn»x2,k, &t4a , &t4b) ; 
terra^OS^aJb (mr * b2 » bm , 1 , &t5a f &t5b) ; 
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> 

t4 = term - 04(u,w, v,bm) ; 

t5 = term_05(u,w,am,bm,v) ; 

return exp ( t l+t2+t3+ (t4a-t4b) *t4/ (R*t ) + (t5a-t5b) * 
t5/(R*t) ) ; 

} 


/* mix_rule.c */ 

/* Fifth (a) in fugacity expression when LMR for "b" */ 

double term„05_a_lmr (bi) 
double bi; 

{ 

return (bi) ; 

} 

/* Fifth (b) in fugacity expression when LMR for "b 3 */ 

double term_05_b_lmr (bm) 
double bm; 

{ 

return bm ; 

} 

/* Fifth (a) in fugacity expression when QMR for "b" */ 

double term_05_a_qmr (bi ,bm, lij ) 
double bi,bm,lij ; 

{ 

return (bi+bm) * (1-lij ) ; 

} 

/* Fifth (b) in fugacity expression when QMR for "b ; */ 

double term_05_b„qmr (bm) 
double bm; 

{ 

return ( 2.0*bm ); 

} 

/* Fourth (b) in fugacity expression when QMR for "a" */ 

double term ...04 _b_qmr ( am) 
double am; 

{ 

return( 2.0*am ); 

> 

/* Fourth (a) in fugacity expression when QMR for "a" */ 
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double term_04_a_qmr (ai , am,kij ) 
double ai,am,kij; 

{ 

return 2. 0*sqrt (ai*am) *(1 . 0 - kij); 

} 

/* Fourth (a) in fugacity expression when CMR for "a" */ 

double term_04_a_cmr (kij ,x2 , am) 
double kij,x2,am; 

{ 

double xl,tl,t2,t3,t4; 
xl = 1.0 - x2; 

tl= pow(mi , . 3333) *pow(a2*mj , .3333) ; 

t2= pow(mj , . 3333)*pov(al*mi, .3333) ; 

t3= sq(xl*pow(mi , .3333) + x2*pow(mj , . 3333) ) ; 

t4= pow(sqrt (al*a2) , .3333) * (1-kij) * (1 - 2*x2) ; 

return( 3*pow(am, . 6667)* ( ((tl-t2)/t3) +t4) ); 

> 

/* Fourth (b) in fugacity expression when CMR for "a" */ 

double term_04„b_cmr (kij ,x2,am) 
double kij,x2,am; 

{ 

double xl,tl,t2,t3,t4; 
xl = 1.0 - x2; 

tl=pow(mi, . 3333) *pow(a2*mj , .3333) ; 
t2=pow(mj , . 3333) *pow(al*mi , .3333) ; 
t3=sq(xl*pow(mi, . 3333) +x2*pow(mj , .3333)) ; 
t4=pow(sqrt (al*a2) , .3333)*(l-kij)*(l-2*(x2*x2+xl* 
xl)); 

return ( 3*pow(am, . 6667) * ( ( (xl-x2) * (t l-t2) /t3) + 
t4) ); 

} 


/* prdeql.c */ 


/*.. Predict Equilibria . .*/ 
extern double avgerror (double) ; 

double predict _eql(ixit optjrnr, int opt_eos, double x[] [MAX J/ARIABLE] ) 
{ 

int i » j ; 

double y2„guess, dif ; 

double guess* 1.0,y2[N] ,diff 00 , sum=0 . 0 ; 

double v^by^SW , am » bm » f ug2_SW » OF , k , 1 > f ug2 JPR ; 

double Q,D f dJ),dJ5 , gamma., i»Ae ; 

double T^pred, P^expCN] ; 

char Pred„Fydat [20] , *aql * "eql.out" ; 



A Program Listing of the Rigorous Model 


86 


FILE *pred_fp,*eql_fp ; 


/★retry :fprintf (stderr, "\nEnter *T-P-y* Filename : ") ; 

scant ('"/,s" ,Pred_Pydat) ;*/ 
if ( (pred_f p=f open(Pydat , "rt") ) == MULL ) 

{ 

fprintf (stderr, "f open failure : 7,s\n" ,Pydat) ; 

> 

f scant (pred_fp , "‘/.lf'/.d" ,&T_pred,£N_exp) ; 
for(i=0;i< N_exp;i++) 

{ 

f scant (pred_fp, r, 7.1f'/,lf " ,&P_exn [i] ,&y2_exp [i] ) ; 

> 

f close (pred_fp) ; 

fprintf (sol.fp, "\n\n** Prediction at 7.21fK **” , 

T_pred) ; 
mr = opt_mr ; 
eos = opt_eos ; 

Select _EoS(eos) ; 

Calc_Constants(eos) ; 

B_1 = bl - al/(R*t) ; 

B_2 = b2 - a2/ (R*t) ; 

k = x [0] [0] ; 
if (no_variables == 1) 

1 = 0.0 ; 
else 

1 = x[0] Cl] ; 

k = 0.0; 

1 = 0 . 0 ; 

for (i=0;i<N_exp;i++) { 
j =0 ; 

y2[i] = y2„exp[i] ; 

do { 

y2„guess= y2[i] ; 
if (mr == 5) { 
eval^Ae(y2.-guess ,P_exp[i] , 

&Ae,&gamma_.i) ; 
eval„D_CLwsmr (y2„guess , 

k , Ae » gamma „ i , &Q , &D , &d w Q , 

fcd.D) ; 

} 

mr^for* a.b ( y 2 .gue ss,k,l,Q,D, &am , 

&bm) ; 

v_by_SW = calc.vCu^w^P^expCi] ,T w pred, 
am, bm, guess) ; 

/* if (!j kk v^by^SW <=» bm ) { 

printf (”\rx bm is greater than v !!!!!”) ; break;} */ 
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fug2_SW = eval_fug_SW(k,l,u,w,P_exp[i] , 
v_by_SW , T_pred , am , bm , D , Q , d_D , d_Q , 
y2_guess) ; 

fug2_PR = eval_f ug_PR_WS (y2_guess , 

P_exp [i] , v_by_SW , T_pred , am , bm , Ae , 
gamma_i,k) ; 

y2[i] = solblty(T_pred,P_exp[i] , 
fug2_SW) ; 

dif = fabs(y2_guess-y2[i]) ; 

> while ( dif > TOLRNS) ; 

> 

for (i=0;i<N_exp;i++) 

diff [i] = (y2[i] - y2_exp[i])*(y2[i] - 
y2_exp [i] ) ; 
sum += diff [i] ; 

> 

OF = avg_ error (sqrt( sum / N_exp) ); 
if ( ! avg_error) { 

printf ("\n Avg. Error is ZERO , Can you believe it !!") ; 
exit (1) ; 

> 

printf ("\n\nAvg. Error of Prediction : ’/.If", OF) ; 

printf ("\n\nCalculated \t Experimental") ; 

for (i=0;i<N_exp;i++) { 

printf ("\n'/.e \t ’/.e" ,y2[i] ,y2_ext>[i]) ; 

> 

f printf (sol_fp , "\n\nCalculated \t Experimental") ; 

f printf (sol_fp,"\n ") ; 

for (i=0; i<N_exp;i++) { 

f printf (sol_fp , "\n'/.e \t */,e" ,y2[i] ,y2_exp[i]) ; 

> 

fprintf (sol_fp, "\n\nAvg. Error of Prediction : */.lf", 

OF) ; 

fprintf (sol_fp , "Xn^****** ********************************* ********\i 1 ") ■ 
return OF ; 

> 


/* rosen.h */ 


/* ■ ■ 

Constant Declarations ..*/ 

#def ine COUNTER 10 
#def ine TOLERENCE .0001 
#def ine MAX.FUNCTION 10 
#def ine MAX.VARIABLE 10 

/*. • 

Global Declarations .,*/ 
double xl'MAX.VARlABLE] [MAX.VARIABLE] ; 

double (#function[MAX_FUNCTIQN] ) 0 ; /* function pointer */ 

int dimension [MAX _V ARI ABLE] , no .test _£unct ion , t imes.f unct ion.calculated ; 
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int no_variables ; 

FILE *out; 

double (*rosen.brock) () , temp_-.-ector [1] [MAX_ V ARI ABLE] ; 
int mr,eos; 


/* rosen.c */ 


#include <stdio.h> 

#include <math.h> 

#include "rosen.h" 

/* extern double opt _paramtrs □ [MAX_ VARIABLE] ;*/ 

void initiate () 

{ 

/*. . 

Total numbers of test function written here ..*/ 
no_test_function=2; 

/*- . 

"func" is the base address of objective the function ..*/ 
function [1] =f unc ; 
dimension [1] =no„variables ; 
function [2] =ycal; 
dimension [2] =no_variables ; 

> 


void sum__matrix(m,l,a,b, c) 
int m,l; 

double a [] [MAX_VARIABLE] , b [] [MAX_VARIABLE] , c [] [MAX_VARIABLE] ; 
{ 

int i,j; 

for(i=0 ; i<m ; ++i) 
f or (j=0 ; j<l ; ++j) 
c[i] [j] =a[i] [j]+b[i] [j]; 
return ; 

} 

void mult ..matrix (m,p,l, a, b, c) 
int m,p,l; 

double a[] [MAX^VARIABLE] ,b [] [MAX_V ARI ABLE] , c [] [INVARIABLE] ; 
{ 

int i , j , k ; 
double sum; 
for(i=0 ; i<m ; ++i) 
for(j=0 ; j <1 ; ++j) 

{ 

sum=0 . 0 ; 
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f or (k=0 ; k<p ; ++k) 
sum+=a[i] [k]*b[k] [j] ; 
c[i] [j]=sum; 

} 

return ; 

> 

void scalar^mult (m, 1 , a , b , constant) 
int m,l; 

double a [] [MAX.VARIABLE] , b G [MAX ^VARIABLE] , constant ; 
{ 

int i , j ; 

for(i=0 ; i<m ; ++i) 
f or ( j=0 ; j<l ; ++j) 
b[i] [j]=a[i] [j] ^constant ; 
return; 

} 


/********************* MAIN PROGRAM ***************/ 

void rosenC) 

{ 

double xnew[l] [MAX J/ARIABLE] , t [1] [MAX_VARIABLE] ; 
double xoldCl] [MAX J/ARIABLE] , c [1] [MAX_VARIABLE] , 
b [1] [MAX.. VARIABLE] ; 
double lamda [MAX_VARIABLE] , tempi ; 
double norm, sum; 
double step_size=0.01,fl,f2; 

double D [MAX^VARIABLE] [MAXJ/ARIABLE] , A [MAX_VARIABLE] [MAX.VARIABLE] ; 

double D^norm,d[MAX_VARIABLE] , change _f lag [MAXJ VARIABLE] ; 

struct 

{ 

double proj [MAX_VARIABLE] ; 

}s [MAX_ VARIABLE] ; 
int stage=0,l,p,n; 

int success [MAX_VARIABLE] failure [MAX .VARIABLE] , 
change_stage_f lag ; 

int min_f lag, n_f unction, j , point , sear cli_dir ; 

if ( (out=f open ("rosen. out" , "w") )— NULL) 

{ 

printf ("\7TJNABLE TO OPEN OUTPUT FILE\n M ) ; 
exit (1) ; 

} 

initiate () ; 

min_flag=l; 
if (min_.flag== , 20 
min_flag=0; 
reenter: 

n_function = 1 ; 

if ((n„function>nojtest_f unction) I I (n_function<= 
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o)) 

{ 

if (no„test_function==0) 

printf("No test Function available \n"); 

exit (1) ; 

goto reenter; 

> 

rosenbrock=f unction [n_f unction] ; /* pointer assignment */ 

n=dimens ion [n__f unction] ; 
for(i=0 ; i<n ; ++i) 

{ 

x[i] [0] =0.0 ; 

} 

for(i=0 ; i<n ; ++i) 
for(j=0 ; j<n ; ++j) 

{ 

if (i==j) 

s[i] .proj [j] =1 . 0 ; 
else 

s [i] . proj [j] =0 ; 

> /* selection of initial direction as principal cordinates dir */ 

step„size = 0.001; 
stage=0; 
point=0 ; 

times^f unction_calculated=0 ; 
f or(j=0 ; j<n ; ++j) 
xnew [0] [j]=x[j] [0] ; 
do 
{ 

+n-stage ; 

for ( j=0 ; j <n ; ++j) { 

t[0].[j]=x[j] [point] ; /* Point may be 0 or anything */ 

} 

for(i=0 ; i<n ; ++i) 
lamda[i] =step_size ; 
search_dir=0; 
for(i=0 ; i<n ; ++i) 

{ 

change_f lag [i] =0 ; 
success [i] =0; 
failure [i]=0; 
d[i]=0; 

} 

do 

{ 

f or (j=0 ; j<n ; ++j) 

{ 

c[0] [j]=s [search_dir] .proj [j] ; 
xold[0] [j]=xnew[0] [j] ; 

> 

scalar _mult ( 1 , n , c , b , lamda [sear ch_dir] ) ; 
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sumjnatrix ( 1 , n , xold , b , xnew) ; 
f l=(*rosenbrock) (xold ) ; 
f 2=(*rosenbrock) (xnew) ; 
if ( ! min., flag) 

{ 

f l=-f 1 ; 

f 2=-f2; 

> 

if(f2>fl) /* Failure occur */ 

{ 

++f ailure[search__dir] ; 
for(j=0 ; j <n ; ++j) 
xnew[0] Cj] =xold[0] [j] ; 
lamda [sear ch_dir] *=“0 . 5 ; 

> 

else 

{ 

++success [search^dir] ; 

++point ; 

d [search_dir] +=lamda [search_dir] ; 
for (j=0 ; j <n ; ++j) 
x[j] [point] =xnew[0] [j] ; 
lamda [sear ch^dir] *=3 . 0 ; 

} 

++search_dir ; 

if (search_dir==n) 

search_dir~=n; 

if (point==MAX_VARIABLE-l) 

point “=(MAX_VARIABLE“1) ; 

for(i=0 ; i<n ; ++i) 

if (success [i] !=0 && failure [i] ! = 

0) 

change_f lag [i] = 

l; 

change_stage_f lag=l ; 
for(i=0 ; i<n ; ++i) 
if (change_stage_f lag) 
if (change_f lag [i] == 

0) 

/* code folded from here */ 
change _stage„f lag= 

°; 

/* unfolding */ 

}while (change _.stage_flag==0) ; 
norm=0 ; 

f or ( j=0 ; j<n ; ++j) 

{ 

templ=xnew[0] [j]-t[0] [j] ; 

templ*=templ; 

norra+=templ; 

> 

norm=sqrt (norm) ; 
if (norm<=T0LERENCE) goto skip; 
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for (1=0 ; l<n ; ++1) 

{ 

f or(j=0 ; j<n ; ++j) 

A Cj] [1] =0; 

f or (i=l ; i<n ; ++i) 

{ 

for(p=0 ; p<n ; ++p) 
c [0] [p] =s [i] . pro j [p] ; 
scalar _mult(l,n,c,b,d[i] ) ; 
f or (j=0 ; j<n ; ++j) 

A Cj] [1] +=b [0] [j]; 

> 

> 


for(i=0 ; i<n ;++i) 

{ 

f or(j=0 ; j<n ; ++j) 
temp_vector CO] Cj] —0 ; 
for (1=0 ; l<i ; ++1) 

{ 

sum=0 . 0 ; 

f or (j=0 ; j<n ;++j) 
stun+=s Cl] .proj Cj] * 

ACj] Ci] ; 

for(p=0 ; p<n ; ++p) 
c CO] Cp] =s Cl] • proj Cp] ; 
scalar _mult(l,n,c,b, stun) ; 
f or ( j=0 ; j <n ; ++j) 
temp_vector CO] Cj]+= 
bCO] Cj] ; 

> 

D_norm=0 ; 

for(p=0 ; p<n ; ++p) 

{ 

DCi] Cp]=ACp] Ci]-temp_vector CO] Cp] ; 
D_norm+=DCi] Cp] *DCi] Cp] ; 

> 

D_norm=sqrt(D_norm) ; 
for(p=0 ; p<n ; ++p) 
s Ci] • proj Cp] =D Ci] Cp] /D_norm ; 

> 

skip : 

printf ("") ; 

}while(norm>=TOLERENCE) ; 


for(i=0 ; i<n ; ++i) 
temp_vector CO] Ci]=xnewC0] Ci] ; 
> 



Appendix B 


Program Listing of the Simplified 
Model 


/* smain.c */ 

#include <stdio.h> 

#include <strings.h> 

#include <string.h> 

#include <math.h> 

#include "rosen.h" 

#include "smply.c" 

#include "rosen.c" 

#include "sol. input . c" 

#include "seql.c" 

extern double avg_error (double) ; 
char *sol_fn = "sol. out" ; 

void main() 

{ 

int i, k = 0, mr.max = 1 ; 
int opt.mr, opt_eos ; 
double dif , pred.value *, 
extern int no.variables ; 
int variable; 
extern char *dat; 

extern double (*rosenbrock) () , temp.vector [] [MAX.VARIABLE] ; 

double Avg.Err, *pred.y2 = vector (6) ; 

double (*fp.prdeql) () ; 

char *opt .rar.name, *opt.eos.name ; 

char *mr.name, *mod.fn = "mod. out" ; 

char *dat.fn = "datanames" ; 

char *prsnt„mr[] = { 

"QMR.&.QMR", 

"CMR.&.QMR", 

"QMR.&.LMR", 
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"CMR_&_LMR" , 
"Wong_&_Sand" 


FILE *mod_fp ; 


if ( (sol_fp = fopen(sol_fn, "at") 

) 

== NULL ) 

{ 

fprintf (stderr , "Error opening file: '/,s\n", 
sol_fn) ; 
exit(l) ; 

} 


if ( (mod_fp = fopen(mod_fn, "at") 

) 

== NULL ) 

fprintf (stderr , "Error opening file: 7,s\n" , 
mod_fn) ; 
exit(l) ; 

> 


if ( (dat_fp = f open(dat_fn, "rt") 

) 

== NULL ) 

fprintf (stderr, "f open failure : 7,s\n" , 
dat_fn) ; 
exit (1) ; 

> 


for (k = 1; k <= 7; k++) 

{ 

/ *tpy _ input () ;*/ 
sy st em_ input () ; 

fprintf (sol_fp, "\n\t*** SYSTEM ***\n") ; 
fprintf (sol_fp, "\n*y.s*\t*‘/.s*\t*'/.s*\n" , 
fn2, fnl, 

Pydat) ; 

fprintf (mod_f p , " Xn'/.sXt’/.sYt'/.sYt " , f n2 , 
fnl , Pydat) ; 

fprintf (sol_fp, " \nTemperature : */,lf\tNo. of Data points : */.d\n" , 
t , N_exp) ; 

fprintf (mod_fp, "'/,.21f\t'/.d" , t, N_exp) ; 

/ *re : */ 
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/★printf ("\nEnter the no of Variables 
scanf ( n 7«d" , too., variables) ;*/ 
no .variables = 2; 

if (no ..variables) { 

f printf (sol_fp, "\n\nThe no. of OPTIMIZATION parameters : 7.d" , 
no_variables) ; 

for (mr = 1; mr <= mr_max; mr++ 

) { 

mr_name = prsnt_mr [ — mr] ; 
mr++ ; 

fprintf (sol_fp, "\n\nMixing Rule for V'aV & V'bV’ is 7,s M , 
mr^name) ; 

for (eos = 2; eos <= 4; 
eos++) 1 

Select „EoS(eos) ; 

Calc^Constants(eos) ; 
rosenQ ; 

Avg_Err = ( ( (*rosenbrock) (temp. vector) ) ) ; 


fprintf (sol.fp, 

M \n7«s ", 

dat) ; 

fprintf (sol.fp, 

"\tkij = 7. .41f \tlij = 7, . 41f \t Avg^Err = %.21f", 
temp.vector [0] [0] , 
temp.vector [0] [1] , 

Avg.Err) ; 
fprintf (mod.fp, 

" \t7. . 41f \t7. . 41f \t7. . 21f " , 
temp ..vector [0] [0] , 
temp.vector [0] [1] , 

Avg^Err * 100.0); 

printf ("\n**AVG. ERROR : 7.1f \n M , 

( ( (*rosenbrock) (temp ..vector) ) ) 

* 100 . 0 ); 


} 

> 

/* goto re;*/ 

> 

} 


/*eql(temp.vector) ; */ 

> 


/* smply.c */ 
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#include "solblty .h" 

#include "mixjrule . c M 

/* Calculation of a,b */ 

void calc_const (eos, tc, pc, accentric, t, omega_a, omega.b, 
a, b) 
int eos; 

double tc, pc, accentric, t, omega_a, omega b, *a, *b; 

{ 

double tr; 
double alpha, beta; 
tr = t / tc; 
switch (eos) { 
case 1: 

*a = omega_a * R * R * tc * tc / pc ; 

*b = omega_b * R * tc / pc ; 
break ; 


case 2: 

*a = omega_a * R * R * pow(tc, 2.0) / pc 

i 

beta = 0.480 + 1.574 * accentric - 0.176 
★accentric * accentric ; 
alpha = 1.0 + beta * (1.0 - sqrt(tr)) ; 
alpha *= alpha; 

*b = omega Jb * R * tc / pc ; 
break; 

case 3: 

★a = omega_a * R * R * tc * tc / pc ; 
beta = 0.37464 + 1.54226 * accentric - 
0.26992 * accentric * accentric ; 
alpha = 1.0 + beta * (1.0 - sqrt(tr)) ; 
alpha *= alpha; 

★a *= alpha; 

*b = omega_b * R * tc / pc ; 
break; 

case 4: 

★a = omega„a * R * R * pow(tc, 2.5) / (pc 
★sqrt(t)) ; 

*b = omega_b * R * tc / pc ; 

break; 

default : 

break; 

> 

} 


/* Schmitt-Wenzel EOS in ’z’ form */ 
double f (A, B, z, u, w) 
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int u, w; 
double A, B, z; 

{ 

return z * z * z + (u * B - B - 1.0) * z * z + 

(w*B*B-u*B*B-u*B + A) * z - (A * 
B+w*B*B+w*B*B*3); 

> 


/* Differential of S-W-Eos */ 

double df(A, B, z, u, w) 
int u, w; 
double A, B, z; 

{ 

return 3.0*z*z+2.0* (u * B - B - 1.0) * 
z+(w*B*B-u*B*B-u*B+A); 

} 


/* Calculation "vap-z" & vap-vol. */ 


double calc_z(u, w, p, t, a, b, guess) 
int u, w; 

double p, t, a, b, guess; 

{ 

double z, diff, dfz, fz, v; 
double A, B; 
int i = 0; 

B = b*p / (R * t) ; 

A = a * p / (R * R * t * t) ; 
z = guess; 
do 
{ 

dfz = df (A, B, z, u, w) ; 

/♦if (dfzCTOL) 

{ 

fprintf (stderr , "Denom < TQL in newton raphson for 
exit (3) ; 

}♦/ 

fz = f (A, B, z, u, w) ; 
guess = z - fz / dfz; 
diff = fabs(z - guess); 
z = guess; 

} while (i++ < ITER kk diff > TOL) ; 
if (diff > TOL) /♦not converged*/ 
return - 100; 

/♦printf ("\nz* # /.lf " ,z) ;*/ 
return z; 

} 
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double fug_PR_infnty (double z, double al, double bt, double 

A, 

double B) 

{ 

double tl, t2, t3, t4, D, rl, r2 ; 

D = u*u-4*w; 
rl = (-u + sqrt(D)) / 2 ; 
r2 = (-u - sqrt(D)) / 2 ; 
tl = bt * (z - 1.0) ; 
t2 = log(z - B) ; 

t3 = A * (2 * al - bt) / (sqrt(D) * B) ; 
t4 = log ( (z - r2 * B) / (z - rl * B) ) ; 

return exp( tl - t2 - t3 * t4 ; ; 

} 


double solblty(t, p, fug2) 
double t, p, fug2 ; 

{ 

double pf ; 

pf = exp(Vs * (p - p_vap) / (R * t)); 

return p_vap * pf / (p * f ug2) ; 

> 

double func(x) 

double x[] [MAX„V ARI ABLE] ; 

{ 

int i, j, stat ; 

double guess = 1.0, y2_cal[N] , diff[N], sum =0.0 

I 

double z_PR, fug2_PR_inf , OF, al, bt , A , B ; 
al = x[0] [0] ; 
bt = x[0] [1] ; 

for (i = 0; i < N_exp; i++) { 

A - al * P[i] / sq(R * t) ; 
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B = bl * P[i] / (R * t) ; 

z_PR = calc_z(u, w, P[i], t, al, bl, guess); 

fug2JPR„inf = fug_PR_infnty (z„?R, al, bt, 

A, B) ; 

y2_cal[i] = solblty(t, P[i], fug2_PR_inf ) ; 

} 

for (i = 0; i < N__exp; i++) { 
dif f [i] = sq( (y2_cal [i] - y2_exp[i]) / 
y2_exp [i] ); 
sum += dif f [i] ; 

} 

OF = sqrt(sum / N_exp) ; 

return OF ; 

> 


double avg_error (OF) 
double OF ; 

{ 

int i; 

double sum = 0.0; 

for (i = 0; i < N_exp; i++) { 

sum += y2_exp[i] ; 

} 

return OF * N_exp * 100.0 / sum ; 

} 


void Select __EoS(eos) 

{ 

switch, (eos) { 

case 1: 

dat = fn[0]; 

break; 

case 2: 

dat = fn[l] ; 

break; 

case 3: 

dat = fn[2] ; 

break; 

case 4: 

dat = fxi[33; 

break; 

default; 
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break ; 

> 

if ( (fp = fopen(dat, "rt")) = JJULL ) { 
fprintf (stderr, "fopen failure : 'Mn", 
dat) ; 
exit (1) ; 

> 

fscanf(fp, "•/.dXdy.lf/.lf &u, kv, &omega_a, &omega_b) ; 
f close(fp) ; 

> 


void Calc_Constants (eos) 

{ 

calc_.con.st (eos , tel, pci, accentricl, t, omega_a, 
omega Jd, &al , &bl) ; 

calc_const (eos , tc2, pc2, accenrric2, t, omega_a, 
omega_b, &a2, &b2) ; 

> 


/* prdeql.c */ 


/*.. Predict Equilibria ..*/ 
extern double avgerr or (double) ; 

void predict_eql(int opt_mr, int opt_eos, double x[] [MAX_VARIABLE] ) 
int i , j , stat ; 

double guess = 1.0, y2_cal[N] , dif f [N] , sum =0.0 
» 

double z_PR, fug2_PR_inf, OF, al, bt , A , B ; 

double T_pred, P_exp[N] ; 

char Pred_Pydat [20] , *eql = "eql.out"; 

FILE * pred_£p, *eql_fp ; 


/♦retry : fprintf (stderr, "\nEnter *T-P-y* Filename : "); 

scanf C"’/.s" ,Pred_Pydat) ;*/ 
if ( (pred.fp = fopen(Pydat, "rt")) == NULL ) { 
fprintf (stderr, "fopen failure : 7.s\n" , 

Pydat) ; 

> 

f scanf (pred_fp, "7.1f'/.d", &T_pred, &N_exp) ; 
for (i - 0; i < N_exp; i++) { 

fscanf (pred_fp, "V.lf/.lf", &P_exp[i] , &y2_exp[i]) ; 
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> 

f close (pred_fp) ; 

fprintf (sol.fp, "\n\n** Prediction at 7,.21fK **", 
T^pred ) ; 

P-vap = vap_pr (T_pred) ; 
mr = opt_mr ; 
eos = opt_eos ; 

Select_EoS (eos) ; 

Calc^Constants (eos) ; 


al = x[0] [0] ; 

bt = x[0] [1] ; 


for (i = 0; i < N_exp; i++) { 

A = al * P[i] / sq(R * t) ; 

B = bl * P[i] / (.R * t) ; 

z_PR = calc_z(u, w, P[i], t, al, bl, guess); 

fug2_PR_inf = fug_PR_infnty(z_?R, al, bt, 

A, B) ; 

y2__cal [i] = solblty(t, P [i] , fue:2_PR_inf ) ; 

> 

for (i = 0; i < N.exp; i++) { 
diff[i] = sq( (y2_cal[i] - y2_exp[i]) / 
y2_exp [i] ); 
sum += diff[i] ; 

} 

OF = sqrt(sum / N_exp) ; 


printf ("\n\nAvg. Error of Prediction : 7, If”, OF) ; 

printf (”\n\nCalculated \t Experimental”) ; 

for (i = 0; i < N_exp; i++) { 

printf (”\n7.e \t 7.e” , y2_cal[i] , y2_exp[i]) ; 

} 

f printf (sol^fp, ”\n\nCalculated \t Experimental”) ; ' 

f printf (sol^fp, ”\n ”) ; 

for (i = 0; i < N„exp; i++) { 
f printf (sol^fp, ”\n7.e \t 7.e”, y2_cal[i], 
y2,.exp[i]) ; 

} 

f printf (sol^fp, ”\n\nAvg. Error of Prediction : 7, If”, 

OF) ; 

f printf (sol^fp , ” \n* ********************************************** *\b- m ) ; 

> 
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/* Sol^input.c */ 

double vap„pr(t) 
double t; 

{ 

if (!strcmp(fn2, "benzoic")) 
return exp(pA2 - pB2 / (t + pC2)) ; 

else if ( ! strcmp(fn2, "anthrcn")) 
return (0.0013332) * exo( (pA2 - pB2 
/ 

(t + pC2) ) ) ; 

else if ( ( ! strcmp(fn2, "acridn")) || ( ! strcmp(fn2, 
"naptqn"))) 

return 1.0e-5 * pow(10.0, (pA2 - pB2 / 
(t+pC2))); 


else 

return pow(10.0, (pA2 - p32 / (t + pC2))); 

} 


/* 

void tpy ..input () 

{ 

fprintf (stderr , "\nEnter *T-P-y* Filename : ") ; 
scanf ("Xs" ,Pydat) ; 

if ( (f d=f open(Pydat , "rt ") ) — NULL ) 

{ 

fprintf (stderr , "f open failure : # /,s\n" ,Pydat) ; 
exit(l) ; 

} 

f scanf (fd,"V.lf */.d y,p_vap M ,&t ,&N„exp,&p_vap) ; 
for(i=0;i< N_exp;i++) 

f scanf (f d, "V.lf'/.lf " ,&P [i] ,&y2_exp[i] ) ; 

} 

f close(fd) ; 

} 

*/ 

void system.*, input () 

{ 

int opt ; 


/♦fprintf (stderr, u \nEnter *SQLVENT* Filename: M ) ; 
scanf (" # /.s'\fnl) ;*/ 



B Program Listing of the Simplified Model 


103 


f scanf (dat_fp, ’7.s 7. s 7 # s\n\ fn2, fnl, Pydat) ; 

if ( (fq = fopen(fnl, "rt")) = NULL ) { 
fprintf (stderr, "fopen failure : 7.s\n", 
fnl) ; 
exit (1) ; 

> 

f scanf (fq, "7,lf 7.1f 7. If 7«lf " , &mi, &tcl, &pcl, &accentricl) ; 
f close(f q) ; 

/ *fprintf (stderr , "\nEnter *T-?-y* Filename : ") ; 
scanf ( 11 y.s ,f , Pydat) ;*/ 

if ( (fd = fopen (Pydat , "rt")) == NULL ) { 
fprintf (stderr, "f open failure : 7s\n" , 

Pydat) ; 
exit (1 ) ; 

} 

f scanf (fd, "7,lf 7.d 7, If", &t, &N_exp, &p_vap) ; 
for (i = 0; i < N_exp; i++) { 
f scanf (fd, "7.1f7.1f", &P[i] , &v2_exp[i]); 

} 

f close(fd) ; 

/*fprintf (stderr , "\n£nter *S0LUTE* Filename : ") ; 
scanf ("7,s" ,fn2) ; */ 

if ( (fr = fopen(fn2, "rt")) == NULL ) { 
fprintf (stderr , "fopen failure : 7,s\n" , 
f n2) ; 
exit (1) ; 

> 

/*printf ("\nDo you have p_sat (1/0):"); 

scanf ("7, d" ,&opt) ; */ 
opt = 1 ; 
if (opt) { 

/*printf ("\nEnter p-sat in bars :"); 

scanf ("7, If H ,&p_vap) ;*/ 
f scanf (fr , "7.1f7,lf7.1f7.1f7.1f ” , taj , &tc2 f 
&pc2, &accentric2, &Vs) ; 

} else { 

f scanf (fr , "7.1f7.1f7.1f 7.1f7.1f7.1f 7.1f7.1f ” , 

&mj , &tc2, &pc2, &accentric2, &Vs, &pA2, 

&pB2, &pC2) ; 
p_vap = vap^pr(t); 

} 

f close(fr) ; 

> 



